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Notation 


The  following  notation  is  used  throughout  this  document.  It  is  primarily  influenced  by  [29] 
as  well  as  having  minor  influence  from  other  authors  listed  in  the  bibliography. 


Vectors 

lower  case  bold  face  letters  such  as  x 

Matrices 

upper  case  bold  face  letters  such  as  P 

Random  Vectors 

lower  case  bold  face  san  serif  type  such  as  x 

Random  Matrices 

upper  case  bold  face  san  serif  type  such  as  P 

Realizations  of 

Random  Vectors  and 
Matrices 

lower  and  upper  case,  respectively,  bold  face 
roman  type  such  as  x  and  P 

Arguments  of 

Densities,  Integrals,  etc. 

Greek  letter  dummy  variables  such  as  ^ 

Superscript  T 

transpose 

Superscript  —1 

inverse 

Superscript  + 

quantity  after  measurement  update  in  the  Kalman  filter 

Superscript  — 

quantity  before  measurement  update  in  the  Kalman  filter 

Raised  “Hat”  (^) 

estimate 

1  •  1 

Matrix  determinant 

Superscript  * 

Quantity  or  variable  associated  with  generalized  residual 
(not  the  complex  conjugate) 
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Abstract 

This  dissertation  develops  a  modification  to  the  standard  Multiple  Model  Adaptive 
Estimator  (MMAE)  which  allows  the  use  of  a  new  “generalized  residual”  in  the  hypothesis 
conditional  probability  calculation.  The  generalized  residual  is  a  linear  combination  of  tra¬ 
ditional  Kalman  filter  residuals  and  “post-fit”  Kalman  filter  residuals  which  are  calculated 
after  measurement  incorporation.  This  modified  MMAE  is  termed  a  Generalized  Residual 
Multiple  Model  Adaptive  Estimator  (GRMMAE).  The  dissertation  provides  a  derivation 
of  the  hypothesis  conditional  probability  formula  which  the  GRMMAE  uses  to  calculate 
probabilities  that  each  elemental  filter  in  the  GRMMAE  contains  the  correct  parameter 
value.  Through  appropriate  choice  of  a  single  scalar  GRMMAE  design  parameter,  the 
GRMMAE  can  be  designed  to  be  equivalent  to  a  traditional  MMAE,  a  post-fit  residual 
modified  MMAE,  or  any  number  of  yet  unused  MMAEs.  The  original  GRMMAE  design 
goal  was  to  choose  the  GRMMAE  design  parameter  which  caused  the  fastest  GRMMAE 
convergence  to  the  correct  hypothesis.  However,  this  dissertation  demonstrates  that  the 
GRMMAE  design  parameter  can  lead  to  /3-dominance,  a  negative  performance  effect  in 
the  GRMMAE.  This  is  a  key  contribution  of  this  research  as  other  researchers  have  previ¬ 
ously  suggested  that  the  use  of  post-fit  residuals  may  be  advantageous  in  certain  MMAE 
applications.  This  dissertation  directly  addresses  the  use  of  post- fit  residuals  by  those  re¬ 
searchers  and  demonstrates  that,  for  their  application,  equivalent  performance  is  achieved 
using  a  traditional  MMAE. 
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Generalized  Residual  Multiple  Model  Adaptive  Estimation 
of  Parameters  and  States 


I.  Introduction 

1.1  Overview 

It  is  often  necessary  to  estimate  both  a  system  parameter  and  the  system  state.  A 
common  example  of  this  need  is  precise  relative  positioning  using  Global  Positioning  Sys¬ 
tem  (GPS)  carrier-phase  measurements.  These  precise  GPS  positions  are  routinely  used 
in  surveying  and  are  planned  for  the  US  Navy  variant  of  the  Joint  Precision  Approach 
and  Landing  System  (JPALS).  The  next  chapter  demonstrates  that  a  constant  parameter, 
namely  carrier-phase  integer  ambiguities,  must  be  estimated  in  order  to  obtain  carrier- 
phase-based  GPS  positions.  The  parameter  estimate  is  then  used  as  part  of  the  system 
model  to  estimate  the  system  state  variables  (typically  including  position,  velocity,  and 
acceleration  in  a  navigation  system).  Thus,  there  is  a  motivation  to  develop  techniques  for 
estimating  states  in  the  presence  of  an  unknown  system  parameter.  In  general,  these  tech¬ 
niques  use  sensor  measurements  to  estimate  the  system  parameter.  In  turn,  the  parameter 
estimate  is  used  to  calculate  the  state  variable  estimates.  A  generic  block  diagram  of  such 
a  system  is  shown  in  Figure  1.1  where  z  is  the  measurement  vector,  x  is  the  state  vector 
estimate,  and  a  is  the  parameter  vector  estimate. 
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Figure  1.1  Generic  Parameter  and  State  Estimator 
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A  second  common  application  for  the  system  in  Figure  1.1  is  failure  detection  and 
isolation.  For  this  application,  the  system  parameter  is  an  operating  mode  (i.e.,  fully 
functional,  type  1  failure,  type  2  failure,  etc.).  The  generic  parameter  and  state  estimator 
uses  sensor  measurements  to  estimate  the  system  mode  and  to  produce  an  appropriate 
state  variable  estimate. 

By  defining  an  appropriate  system  parameter,  both  of  these  common  estimation  ap¬ 
plications  can  be  addressed  using  the  generic  system  of  Figure  1.1.  Since  this  research  was 
inspired  by  GPS  carrier-phase  integer  ambiguity  resolution,  that  application  is  discussed 
in  some  detail.  However,  this  research  is  sufficiently  general  to  apply  to  the  failure  detec¬ 
tion  and  isolation  application  as  well  as  other  appropriate  applications  where  the  system 
state  is  estimated  in  the  presence  of  an  unknown  system  parameter.  Chapters  II  and  III 
explicitly  state  all  mathematical  assumptions  an  application  must  meet  for  this  research 
to  apply. 

The  Multiple  Model  Adaptive  Estimator  (MMAE)  [26]  of  Figure  1.2  is  a  specific  type 
of  parameter  and  state  estimator.  The  MMAE  consists  of  a  bank  of  parallel  Kalman  filters, 
termed  “elemental  filters.”  Each  elemental  filter  is  based  on  a  hypothesized  parameter 
vector,  ai,  a2,  ...  or,  slk,  where  aj  is  the  constant  (or  slowly- varying  with  respect  to 
the  states)  parameter  vector  for  the  elemental  filter,  and  there  are  K  elemental  filters 
in  the  MMAE.  The  elemental  filters  act  upon  a  measurement  vector  z  to  produce  a  state 
estimate  Xj,  a  residual  vector  rj,  and  a  residual  covariance  matrix,  Aj.  The  MMAE  is 
derived  in  detail  in  Chapter  II. 

The  primary  motivation  for  using  a  multiple  filter  technique  like  the  MMAE  rather 
than  a  single  filter  technique  is  responsiveness.  In  many  cases,  problems  that  are  addressed 
using  the  MMAE  could  also  be  addressed  with  a  single  filter  estimation  technique  such  as 
a  Kalman  filter.  However,  by  comparing  measurement  residuals,  r,  from  several  elemental 
filters,  the  MMAE  is  generally  able  to  determine  the  correct  parameter  value  more  quickly. 
Additionally,  the  MMAE  will  generally  respond  to  changes  in  the  parameter,  such  as  a  sys¬ 
tem  failure,  more  quickly  than  a  single  filter  solution  [46].  The  additional  filters  required 
in  the  MMAE  do  increase  the  computational  load  when  compared  to  a  single  filter  solu¬ 
tion.  However,  computational  power  continues  to  increase,  and  computationally  efficient 
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Figure  1.2  General  MMAE  Block  Diagram 


algorithms  such  as  “moving-bank”  and  “filter  spawning”  MMAE’s  have  been  developed  to 
address  this  disadvantage  [10,  11,  12,  13,  14,  23,  24,  28,  41,  42,  46].  This  research  assumes 
that  sufficient  computational  power  exists,  so  methods  to  increase  computational  efficiency 
are  not  emphasized. 

Recently,  researchers  have  suggested  modifying  the  MMAE  to  calculate  hypothesis 
conditional  probabilities  from  “post-fit”  residuals  [17,  18].  The  traditional  MMAE  uses  a 
nonlinear  function  of  measurement  residuals  from  each  elemental  filter  to  calculate  prob¬ 
abilities  that  the  filter  contains  the  correct  hypothesis.  These  researchers  have  suggested 
that  better  performance  may  be  achieved  by  calculating  hypothesis  conditional  probabil¬ 
ities  using  a  nonlinear  function  of  residuals  calculated  after  measurement  update.  Such 
residuals  are  termed  “post-fit”  residuals.  Additionally,  [47]  has  suggested  using  post-fit 
residuals  from  multiple  Kalman  filters  for  hypothesis  testing.  However,  it  is  unclear  from 
[47]  if  a  full  MMAE  was  used  or  if  post-fit  residual  monitoring  was  used.  Regardless  of 
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what  technique  [47]  used,  the  method  is  at  least  closely  related  to  using  post-fit  residuals 
in  an  MMAE. 


This  dissertation  was  inspired  by  [17,  18],  which  provided  no  formal  derivation  or 
mathematical  characterization  of  the  effect  that  post-fit  residual  processing  has  on  the 
MMAE.  A  research  contribution  of  this  dissertation  is  to  fill  the  gap  left  in  [17,  18].  To  do 
this,  a  mathematical  generalization  to  the  MMAE  is  developed  to  allow  the  designer  easily 
to  incorporate  linear  combinations  of  Kalman  filter  residuals  (calculated  before  measure¬ 
ment  update)  and  post-fit  residuals  (calculated  after  measurement  update).  This  linear 
combination  is  termed  a  “generalized  residual.”  The  MMAE  acting  upon  this  generalized 
residual  is  introduced  here  as  the  “Generalized  Residual  Multiple  Model  Adaptive  Estima¬ 
tor”  (GRMMAE).  This  research  provides  a  derivation  of  the  hypothesis  conditional  prob¬ 
ability  calculation  based  upon  the  generalized  residual.  Additionally,  generalized  residual 
properties  for  both  correct  and  incorrect  hypotheses  are  derived.  The  effect  of  using  the 
generalized  residual  on  MMAE  performance  is  demonstrated.  The  fact  that  the  general¬ 
ized  residual  leads  to  the  /3-dominance  effect  [2,  3,  32],  a  negative  performance  attribute, 
is  demonstrated.  Additionally,  for  the  specific  design  and  application  used  by  [17,  18,  21], 
this  dissertation  demonstrates  that  the  /3-dominance  effect  is  cancelled  and  proves  that  the 
suggestion  in  [17,  18]  that  post-fit  residuals  improve  performance  for  GPS  carrier-phase 
integer  ambiguity  resolution  is  not  true.  These  are  all  contributions  of  this  research.  Eur- 
ther,  this  dissertation  demonstrates  that,  in  general,  the  proposed  GRMMAE  is  susceptible 
to  the  /3-dominance  effect  and  is  therefore,  not  recommended  for  wide  application.  How¬ 
ever,  as  the  recommendation  for  future  research  suggests,  with  further  development,  the 
GRMMAE  may  be  useful  to  correct  for  /3-dominance  optimally. 

1.2  Dissertation  Summary 

This  dissertation  consists  of  5  chapters.  This  chapter  provides  an  introduction  to 
the  dissertation.  Ghapter  II  gives  necessary  background  for  the  new  research.  Ghapter  III 
derives  the  new  theory  for  this  dissertation,  presenting  a  full  derivation  of  the  proposed 
GRMMAE  and  properties  of  the  generalized  residuals.  The  primary  research  contribu¬ 
tions  of  this  dissertation  are  contained  in  Ghapter  III.  Ghapter  IV  presents  simulation  and 
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analysis  results  for  the  GRMMAE,  demonstrating  that  the  modified  MMAE  from  [17,  18] 
(a  specific  example  of  the  GRMMAE)  provides  equivalent  performance  to  the  traditional 
MMAE.  Additionally,  it  demonstrates  that,  in  general,  the  proposed  GRMMAE  is  sus¬ 
ceptible  to  the  /3-dominance  effect.  Simulation  results  in  Ghapter  IV  validate  theoretical 
contributions  presented  in  Ghapter  III.  Ghapter  V  gives  conclusions  and  recommenda¬ 
tions  based  on  this  research.  Additionally,  the  research  contributions  of  this  dissertation 
are  enumerated  in  Ghapter  V. 
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II.  Background 


2.1  Kalman  Filter 

The  basic  building  block  in  the  Multiple  Model  Adaptive  Estimator  (MMAE)  is  the 
Kalman  filter.  Therefore,  the  Kalman  filter  is  derived  here.  This  section  presents  the  linear 
Kalman  filter  used  for  systems  which  can  be  described  by  linear  dynamics  and  measurement 
models.  Section  2.2  describes  the  extended  Kalman  filter,  which  is  used  when  a  system 
has  nonlinearities  in  the  dynamics  and/or  measurement  model. 

Consider  a  system  which  can  be  modeled  by  the  linear  stochastic  difference  equation 


x(tj+i)  =  +  Bd{ti)u{ti)  +  Wd{ti) 

z{ti)  =U{ti)x{ti)  +  v(ti)  (2.1) 

In  the  first  line  of  Equation  (2.1),  x  (t*)  is  the  state  vector,  $  (tj+i,  tj)  is  the  state  transition 
matrix  from  time  ti  to  ti+i,  Bd{ti)  is  the  input  matrix,  u(tj)  is  the  control  vector,  and 
Wrf  (ti)  is  the  discrete-time  dynamics  driving  noise  vector.  This  represents  the  system 
dynamics  model.  The  second  line  of  Equation  (2.1)  gives  the  system  measurement  model, 
where  z  (ti)  is  the  measurement  vector,  H  (ti)  is  the  measurement  matrix,  and  v  (ti)  is  the 
measurement  noise  vector.  Note  that  Equation  (2.1)  may  be  an  equivalent  discrete-time 
description  of  a  continuous-time  system  model  [30].  The  form  of  Equation  (2.1)  assumes 
measurements  z  (ti)  are  available  only  at  discrete-times. 

Noise  vectors  Wd  (ti)  and  v  (ti)  are  assumed  independent,  zero-mean,  white,  Gaussian 
processes  with  covariance  kernels 


E  {\Md  (ti)  wj  (tj)}  =  Qd  (ti)  kj 

E{v(ti)  v^(tj)}  =  B(ti)5ij  (2.2) 

where  the  discrete-time  dynamics  driving  noise  covariance  (ti)  is  assumed  positive  semi- 
definite  for  all  time,  the  measurement  noise  covariance  R  (ti)  is  assumed  positive  definite  for 
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all  time,  and  5ij  is  the  Kronecker  delta  function.  It  is  assumed  the  system  of  Equation  (2.1) 
starts  from  an  initial  state  x(fo),  which  is  a  Gaussian  random  vector  with  mean  and 
covariance  given  by 


E{x(to)}  =  xo 

E{[x(to) -xo]  [x(to) =  Po  (2.3) 

where  xq  is  the  initial  state  vector  estimate,  x  (to)  is  assumed  independent  of  (tj)  and 
V  (tj)  for  all  ti,  and  Pq  is  the  initial  state  covariance  matrix.  Pq  is,  by  definition,  positive 
semi-definite. 

Two  measurement  times  U-i  and  ti  are  used  to  derive  a  two-step  process.  The  first 
step  propagates  the  optimal  state  estimate  and  state  covariance  from  the  time  immediately 
after  measurement  update  at  ti-i  to  the  time  immediately  before  measurement  update 
at  ti-  The  next  step  updates  the  propagated  state  and  covariance  estimates  using  the 
measurement  update  at  time  ti.  Notationally,  t~  is  used  to  denote  a  quantity  at  time  U 
immediately  before  the  measurement  update  is  incorporated,  while  tf  denotes  a  quantity 
immediately  after  the  measurement  update  at  time  ti.  The  derived  propagate  and  update 
algorithm  for  the  general  times  ti-i  and  ti  is  used  recursively  to  estimate  the  state  vector 
and  state  covariance,  beginning  at  the  initial  conditions  of  Equation  (2.3). 

Suppose  the  measurement  z  (U-i)  has  been  incorporated.  The  desire  is  to  define 
the  conditional  probability  density  function  (PDE)  for  the  state  vector  conditioned  on 
the  entire  measurement  history  /x(q_i)|z(ti_i)  (C|Z  (ti-i)  =  and  to  propagate  that 

density  forward  until  the  next  measurement  is  incorporated.  Here,  Z(tj_i)  is  a  growing 
length  measurement  history  vector  consisting  of  all  measurement  vectors  from  to  through 
ti-i,  ^  and  Zi-i  are  dummy  variables  representing  the  state  vector  and  the  measurement 
history,  respectively.  It  is  proven  in  [29]  that /x(t._,^)|z(q_i)  (4|Z  (ti-i)  =  .Zj_i)  is  Gaussian 
under  the  previous  assumptions  concerning  Equations  (2.1)  through  (2.3).  Therefore,  the 
mean  and  covariance  fully  define  the  density.  Eor  convenience,  define 
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(2.4) 


X  {tf-i) 

p  {4-i) 


x(ti_i)  I  =  Zi_i 


x(ti_i)  -  X  (t+ ^)]  [x(ti_i) 


S(tti)]^  I  Z(ti-i) 


Zj-i 


where  x  and  P  are  the  conditional  state  mean  and  covariance  estimates,  con¬ 

ditioned  on  the  realized  measurement  history  through  time  tj-i,  Zj_i,  immediately  after 
measurement  incorporation  at  ti-i.  For  a  Gaussian  density,  the  conditional  mean  x  (^(^1) 
is  the  optimal  state  estimate  following  measurement  update  at  time  ti-i  for  “essentially  any 
criterion  of  optimality”  [29].  Not  only  is  it  the  conditional  mean,  it  is  also  the  conditional 
mode  and  the  conditional  median. 

Next,  the  conditional  mean  and  covariance  need  to  be  propagated  to  time  t~ .  To 
simplify  this  discussion,  assume  there  are  no  control  inputs  (i.e.,  u(tj_i)  =  0).  Control 
inputs  only  have  the  effect  of  shifting  the  state  conditional  mean.  They  have  no  effect  on 
the  conditional  covariance.  The  effect  of  control  inputs  is  added  following  this  development. 
The  state  conditional  mean  at  ti,  conditioned  on  the  measurement  history  at  is 


E 


x{ti)  I  Z{ti-i)  =  Zj_i 


=  E 


^  -h  Wrf(ti_i)  \Z{ti-i)  =  Zi_i 


=  ^{ti,ti-i)E  \Z{ti-i)  =  Zi_] 


+  E 


Wd  (tj-l 


E 


x{ti)  \Z{ti-i)  =  Zj_i  =  ^{ti,ti-i)E  x{ti-i)  I  Z(ti_i)  =  Zi_i 


(2.5) 


In  the  third  line  of  Equation  (2.5),  the  unconditional  mean  is  used  because  Wrf(tj)  is 
independent  of  Z(tj_i).  In  the  final  line,  the  assumption  that  w^(tj)  is  zero-mean  is 
invoked. 

Define 
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x(t.  )  =  E  x{ti)  I  =  Zi_i 

[^(*»)  -  I  Z(t*-i)  =  Zi_i  (2.6) 

Using  Equations  (2.5)  and  (2.6),  the  propagated  state  conditional  mean  is 

^{ti)  =  ^(U,U-i)S  (tti)  (2-7) 

The  propagated  state  covariance  estimate  is 

[^(^*)  -  ^(^*“)]  [^(*»)  -  I  Z(U-i)  =  Zi_i 

=  E  [$  (U,U_i)x(U_i)  +  Wrf(U)  -  x(t“)]  [$(U,ti_i)x(U_i)  +  \Md{ti) 

-  x(t“)]^  I  Z(U_i)  =  Zi_i 

=  ^{ti,ti-i)E  x(U_i)x^(U_i)  I  Z(U_i)  =  Zi_i  $'^(ti,U_i) 

-x(t“)E  x'^(U_i)  I  Z(U_i)  =  Zi_i  $'^(ti,U_i) 

-  ^  E  x(U_i)  I  Z(U_i)  =  Zi-i  x^(t“) 

+  x(t~)x^(t~)  +  E  Wrf(U)wJ(ti)|  Z(ti_i)  =  Zi_i 
=  $(U,U_i)E  x(U_i)x^(U_i)  I  Z(U_i)  =  Zi_i  $'^(U,U_i) 

-  X  (t“ )  (tti)  ti-i)  -  $  (U,  U-i)  X  (t+i)  (^*“ )  +  X  [t- )  x^  {t- ) 

+  Qd  (U) 

=  ^{ti,ti-i)E  x(U_i)x^(U_i)  I  Z(U_i)  =  Zi_i  $'^(U,U_i) 

-  $  (ti,  ti-l)  X  X^ 

-  $  (u,  ti-l)  X  (t+ i)  X^  (u,  U-i) 

+  $(U,ti_i)x(t+i)x^(t+i)^'^(U,U-i)  +  Qd{ti) 
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P  <  E  X  (ti-i)  (ti-i)  \Z{ti-i)  =  Zi_i 


+  Qd{ti) 

P(t“)  =  ^  ^  +  Qd{ti) 


(2.8) 


where  the  cross-terms  involving  x  (ti)  with  (U)  were  dropped  because  Wd  (U)  is  assumed 
to  be  independent  of  x(tj)  and  zero-mean.  Similarly,  the  cross-terms  with  and 

Wrf  (ti)  were  dropped  because  (t*)  is  zero-mean. 

If  X  (t~ )  is  used  as  the  optimal  propagated  state  estimate  prior  to  measurement 
update  at  ti,  then  the  error  associated  with  this  estimate  is  [x(tj)  —  x  Since  this 

error  is  zero-mean  [29] ,  P  {t~ )  is  the  estimation  error  covariance  as  well  as  the  propagated 
state  estimate  covariance.  Note  also  that  the  measurement  vector  z  (tj)  does  not  appear  in 
Equation  (2.8).  Therefore,  given  a  covariance  for  the  initial  state  estimate,  the  propagated 
state  estimate  covariance  (and  state  estimate  error  covariance)  can  be  computed  for  all 
time  before  any  measurements  are  taken. 

Now,  consider  updating  the  propagated  state  estimate  with  information  from  the 
measurement  taken  at  ti.  Thus,  the  conditional  density  /x(q)|z(ti)  (^|Z  (ti)  =  tZi)  is  de¬ 
sired.  Repeated  application  of  Baye’s  rule  and  the  definition  of  the  measurement  history 
yield  [29] 


/x(q)|Z(ti)  (^|Z(ti)  =  Zi) 

_  /z(ti)|x(ti),Z(q_i)  (Cil^  (ti)  z  (ti—l)  Zi—i)  (^jZ  (ti— l) 

/z(ti)|Z(ii_i)  (CijZ  (tj—i)  Zi—i) 


Z^-i) 


(2.9) 


where  Ci  is  the  dummy  variable  associated  with  z  (ti). 

The  form  of  Equation  (2.9)  is  motivated  by  the  fact  that  all  three  densities  on  the 
right  hand  side  are  easily  evaluated.  The  density  /x(q)|z(ti_i)  (C|Z(ti-i)  =  Zi-i)  is  the 
density  for  the  propagated  state,  developed  earlier  in  this  section.  It  is  Gaussian  [29],  with 
mean  and  covariance  given  by  Equations  (2.7)  and  (2.8): 
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/x(ti)|Z(ti_i)  (^1^  (^i— l)  ^i—l) 

“  (2^)"/^  IP  (t-)  I V2  ”P  ^  (*r)]^p~’  (ti )  [{ -  g  (ti )]  }  (2.10) 

where  n  is  the  dimension  of  the  state  vector  and  |  •  |  is  the  matrix  determinant. 

Next,  consider  the  first  conditional  density  in  the  numerator  of  Equation  (2.9), 
fz{ti)\x(ti),z{ti_i)  iCMiU)  =  =  Zi_i).  Recall  from  Equation  (2.1)  that  the  mea¬ 

surement  model  is 


z{ti)  =  H(tj)x(ti)  -h  v(ti) 

Since  this  density  is  conditioned  on  the  fact  the  state  has  taken  a  value,  x  (tj)  = 
then  the  random  vector  z  (ti)  is  a  linear  function  of  the  known  state  realization  and  the 
zero-mean,  white,  Gaussian  measurement  noise  with  covariance  R(ti).  Therefore,  this 
conditional  density  is  a  Gaussian  density,  with  mean  H  (tj)  ^  and  covariance  R  (tj): 

/z(ti)|x(ii),Z(ti_i)  (Cil^  (^i)  Cl  Z  (tj_i)  Zi—i) 

where  m  is  the  dimension  of  the  measurement  vector. 

Einally,  consider  /z(q)|z(ti_i)  (Ci|Z  (ti_i)  =  .Zj-i),  from  Equation  (2.9),  which  is 
Gaussian  with  mean  H  (ti)  x  and  covariance  ^{h)  (ii)  P  ih  )  (*»)  +  P 

calculated  from  the  measurement  model  in  Equation  (2.1)  [29]: 

/z(ti)|Z(q_i)  (Cil^  (^i— l)  ^i—l) 

=  "  H(tOX(t-)]^A-i(t-)  [C  -  H(tOX(t-)] 

(2.12) 
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Substituting  Equations  (2.10)  through  (2.12)  into  Equation  (2.9),  and  invoking  the 
matrix  inversion  lemma,  yields  the  equivalent  Gaussian  density  [29] 


/x(ti)|Z(ii)  (^1^  (^i)  ~  ^i) 

=  (2„)~/2|p  |i/'2  ~  ^  (2-13) 

where  the  updated  state  vector  mean  and  covariance  immediately  after  measurement  in¬ 
corporation  are 

X(t+)  =X(t-)  +  P{t-)U^{U)[U{U)P{t-)U^{U)  +  R(t.)]“'[z,  -  H(t,)S(t-)] 
P(t+)  =P{t-)  -  P(t-)H^(t2)[H(t,)P(t-)H^(t,)  +  K{U)]-^U{U)P{t-) 

(2.14) 


and  Zj  is  the  realized  measurement  at  time  ti. 

The  Kalman  filter  propagate  and  update  algorithm  has  now  been  developed  assuming 
no  control  inputs  are  present.  To  complete  the  development,  the  effect  of  control  inputs 
must  be  included.  As  mentioned  previously,  these  deterministic  inputs  only  affect  the  mean 
of  the  propagated  state  estimate.  When  control  inputs  are  included,  the  propagated  state 
vector  estimate  becomes  [29] 


x(t.  )  =  $  (ti,tj_i)x  +  Brf(ti)u(ti)  (2.15) 

where  the  control  input  terms,  from  Equation  (2.1),  are  the  equivalent  discrete-time  control 
model  calculated  from  the  continuous-time  model  using 

U 

Bditi)u{ti)  =  j  $(ti,r)B(T)u(r)(ir  (2.16) 

U-i 
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The  Kalman  filter  is  summarized  by  the  propagation  and  update  algorithm  shown  in 
Table  2.1.  In  Table  2.1,  K  (tj)  is  the  Kalman  filter  gain.  All  other  variables  and  equations 
are  previously  defined  in  this  section. 

2.2  Extended  Kalman  Filter 

Frequently,  applications  call  for  use  of  a  Kalman  filter  on  a  nonlinear  system.  Such 
a  system  may  contain  nonlinearities  in  the  dynamics  and/or  the  measurement  model.  One 
example  with  a  nonlinear  measurement  model  is  estimation  of  position  using  the  GPS. 
A  common  solution  for  these  nonlinear  problems  is  the  extended  Kalman  filter.  Like  the 
Kalman  filter  of  Section  2.1,  the  extended  Kalman  filter  can  be  a  basic  building  block  in  the 
MMAE.  Therefore,  the  extended  Kalman  filter  is  derived  here,  paralleling  the  derivation  in 
[30].  It  is  important  to  note,  however,  that  the  MMAE  based  on  extended  Kalman  filters 
is  an  ad  hoc  form.  Unlike  the  MMAE  based  on  linear  system  models  and  linear  Kalman 
filters,  it  is  not  an  optimal  Bayesian  form. 

Consider  a  system  which  can  be  modeled  by  the  continuous-time  dynamics  and 
discrete-time  measurement  models 


X  (t)  =  f  [x  (t) ,  u{t) ,  t]  -h  G  (t)  w  (t) 

z{ti)  =h[x(L),L]  -F  v(L)  (2.17) 


Table  2.1  Kalman  Filter  Algorithm  Summary  [29] 


State  Propagation 

s  (*i  )  = 

^  {ti,ti-i)x{tf_^)  -F  Bd(L)u(L) 

= 

$  (L,L_i)  P  $  (L,L_i)  -F  Qd{ti) 

Measurement  Update 

K  (ti)  = 

P(t-)H^(L)[H(L)P(t-)H^(L)  +  R(L)]  ' 

^{tf)  = 

+  K(L)[zi  -  H(L)X(t“)] 

^{tf)  = 

P(f-)  -  K(L)H(L)P(t-) 
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where 


X  (t)  is  the  state  vector  derivative  at  time  t 

f  [•,  •,  •]  is  the  general,  possibly  nonlinear,  dynamics  model;  a  function  of  x,  u 
and  t 

h  [•,  •]  is  the  general,  possibly  nonlinear,  measurement  model;  a  function  of  x 
and  ti 

For  this  derivation,  it  is  assumed  that  f  [x  (t)  ,  u  (t)  ,  t]  is  Lipschitz  in  x  (t),  continuous  in 
u  (t),  and  piecewise  continuous  in  t.  Additionally,  the  control  input  u  (t)  is  assumed  to  be 
piecewise  continuous  in 

Assume  w  (t)  and  v  (tj)  are  independent,  zero-mean,  white,  Gaussian  noise  processes 
with  covariance  kernels 


-h  r)}  =Q(t)(5(r) 

E  {v  {ti)  {tj)}  ='R{ti)6ij  (2.18) 

where  Q  {t)  is  the  continuous-time  dynamics  driving  noise  strength  and  5  (r)  is  the  Dirac 
delta  function.  Finally,  the  model  is  assumed  to  start  from  initial  state  x(to),  which  is 
a  Gaussian  random  vector  with  mean  xq  and  covariance  Pq.  Note  that  in  the  case  of  a 
nonlinear  dynamics  model  f  [x  {t)  ,  u{t)  ,  t],  the  state  vector  generally  does  not  remain 
Gaussian  beyond  the  initial  time. 

As  a  stepping  stone  to  deriving  the  extended  Kalman  filter,  first  consider  the  lin¬ 
earized  Kalman  filter.  The  linearized  Kalman  filter  assumes  a  nominal  state  trajectory 
x„  {t)  is  defined  for  all  times  of  interest  (defined  as  t  G  T).  Linear  perturbation  theory 
is  then  used  to  calculate  state  trajectory  perturbations  from  the  nominal  trajectory.  This 
theory  is  derived  first,  then  the  extended  Kalman  filter  is  related  to  it. 

Assume  the  nominal  state  trajectory,  subject  to  the  initial  condition  x^  (to)  =  x„q, 
can  be  generated  such  that  it  satisfies  the  deterministic  differential  equation: 

^  These  assumptions  are  necessary  to  ensure  that  a  solution  exists  to  the  deterministic  differential  equa¬ 
tion  X  (t)  =  f  [x  (t)  ,  u  (t) ,  t]  subject  to  the  initial  condition  x(to)  =  xq.  For  more  details,  see  [29] 
Section  2.3. 
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x„  (t)  =  f  [x„  (t) ,  u  (t)  ,  t] 


(2.19) 


To  maintain  the  adequacy  of  linear  perturbation  techniques,  perturbations  of  the  control 
vector  u  (t)  from  a  nominal  control  vector  Ujj  (t)  are  not  considered  here. 

This  nominal  state  process  has  a  set  of  nominal  noise- free  measurements  (tj)  as¬ 
sociated  with  it,  given  by 


Zn  {ti)  =  h  [x„  {ti)  ,  ti]  (2.20) 

In  Equations  (2.19)  and  (2.20),  f  [•,  •,  •]  and  h  [•,  •]  are  the  general  dynamics  and  measure¬ 
ment  models,  respectively,  given  in  Equation  (2.17). 

Consider  the  state  perturbation  from  the  nominal  state.  Using  Equations  (2.17)  and 
(2.19),  the  state  perturbation  from  the  nominal  state  is  given  by 


[x(t)  -  x„(t)]  =  f[x(t),  u{t),t]  -  f[x„(t),  u{t),t]  +  G{t)w{t)  (2.21) 


Note  in  Equation  (2.21)  that  x  (t)  and  x  (t)  are  random  vectors  while  x„  (t)  and  Xn  (t)  are 
deterministic. 

To  develop  the  perturbation  equation,  the  right  hand  side  term  in  Equation  (2.21)  is 
expanded  using  a  Taylor  series  expansion  about  the  nominal  state  vector  Xjj  (t),  yielding 


f[x(t),  u(t),t] 


f  [x„  (t) 


u  (t)  ,  t]  + 


9f[x,  u{t),t] 
dx 


[x  (t) 

X=Xn 


x„(t)]  -h  h.o.t. 

(2.22) 


In  Equation  (2.22),  h.o.t.  represents  higher  order  terms  and  the  derivative  is  assumed  to 
exist.  It  is  assumed  that  the  state  perturbations  from  nominal  are  sufficiently  small  to 
allow  h.o.t.  to  be  neglected.  Substituting  Equation  (2.22)  into  Equation  (2.21)  yields  a 
first-order  linear  approximation  to  the  state  perturbation  from  nominal 
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6x.{t)  =  F[t;x„(t)]  5x(t)  +  G(t)w(t) 


(2.23) 


where  Sx  {t)  is  a  first-order  approximation  of  [x  (t)  —  x„  (t)]  and  F  [t;  x„  (t)]  is  a  matrix 
of  partial  derivatives  of  the  dynamics  model  defined  by 


F[t;x„(t)]  = 


(9f  [x,  u  (t) ,  t] 
3x 


(2.24) 


Thus,  a  linear  approximation  for  the  perturbation  of  state  vector  from  the  nominal  state 
vector  is  defined.  To  calculate  the  total  state  vector  estimate,  the  perturbation  estimate 
(to  be  derived)  is  added  to  the  nominal  state:  x(t)  =  <5x  (t)  -|-  x„(t). 

The  same  approach  used  to  derive  Equation  (2.23)  is  used  to  derive  the  measurement 
perturbation 


5z{ti)  =  H[ti;x„(ti)]  dx(ti)  -h  v(fi) 


(2.25) 


As  before,  Sz  [ti]  is  a  first-order  approximation  to  [z  (tj)  —  (tj)]  and  H  [tj;  x„  (tj)]  is  a 

matrix  of  partial  derivatives  (assumed  to  exist)  of  the  measurement  model  defined  by 


H  \tii  (^i)] 


A 


[x,  tj] 

dx 


x=x„(q) 


(2.26) 


Equations  (2.23)  and  (2.25)  yield  a  linear  model  (albeit  an  approximation)  to  which 
linear  filtering  theory  can  be  applied.  This  form  of  Kalman  filter  is  called  a  linearized 
Kalman  filter  or  perturbation  Kalman  filter.  The  extended  Kalman  filter  is  directly  related 
to  the  linearized  Kalman  filter.  The  difference  between  the  two  lies  in  the  fact  that  the 
extended  Kalman  filter  re-linearizes  about  the  state  estimate  after  each  measurement  up¬ 
date,  rather  than  linearizing  about  a  pre-defined  nominal  state  trajectory.  The  extended 
Kalman  filter  has  the  advantage  over  the  linearized  Kalman  filter  of  being  more  robust  in 
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the  presence  of  large  perturbations  from  the  nominal  state  trajectory.  For  the  extended 
Kalman  filter  to  operate  properly,  though,  any  perturbations  between  two  consecutive 
measurement  updates  must  be  small  enough  to  allow  h.o.t.  from  Equation  (2.22)  to  be 
neglected. 

To  derive  the  extended  Kalman  filter,  assume  a  measurement  update  has  taken  place 
at  time  tj,  producing  state  estimate  x  (tf ) .  Now,  rather  than  maintaining  the  linearization 
about  the  nominal  state  re-linearize  about  the  new  state  estimate.  Let  Xn{t/ti) 

be  the  extended  Kalman  filter  “nominal”  state  trajectory  on  the  interval  Thus, 

^nit/ti)  is  the  solution  to  Equation  (2.19)  subject  to  initial  condition  x(t)'“),  solved  on 
the  interval  [ti,L+i): 


±n{t/ti)  =  f  [x„(t/L)  ,u(ti)  ,f] 

for  (2.27) 

Xn(L/L)  =  x(t+) 

where  f  [•,  •,  •]  is  the  dynamics  model  from  Equation  (2.17)  and  x  {tf )  is  the  state  estimate 
after  measurement  update  at  time  ti. 

To  propagate  the  state  to  the  next  measurement  time  tj+i,  employ  the  re- linearized 
F  [t]  x„  {t/ti)]  where  F  [•  ;  •]  is  defined  in  Equation  (2.24).  Let  Sx  {t/ti)  denote  the  estimate 
of  dx{t)  based  on  measurements  through  time  ti.  Then,  on  the  interval  [ti,ti+i) ,  5x{t/ti) 
is  the  solution  to 

Sx{t/ti)  =  F[t;Xn{t/ti)]Sx{t/ti)  (2.28) 

subject  to  the  initial  condition 

5x{ti/ti)  =  Sx{tf)  =  0  (2.29) 

This  initial  condition  comes  from  the  fact  that  the  best  estimate  of  5x  (ti)  immediately 
after  measurement  incorporation  is  0  due  to  the  re-linearization.  Solving  Equation  (2.28) 
subject  to  Equation  (2.29)  leads  to  Sx(t/ti)  =  0  over  the  interval  [ti,L+i).  Thus,  the 
propagated  state  perturbation  estimate  fe  (^i+i)  is  0. 
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Next,  the  filter  is  updated  with  a  measurement  as  in  Equation  (2.25),  with  z„  (ti+i) 
replaced  by 

Zn  {ti+i/ti)  =  h  [x„  {ti+i/ti) ,  ti+i]  (2.30) 

where  Zjj  (tj+i/tj)  is  the  measurement  prediction  based  on  re-linearization  about  the  most 
recent  state  estimate.  Given  the  measurements  of  (2.25)  and  (2.30),  the  measurement 
update  of  the  state  perturbation  estimate  is 


<5x  {t 


+  K(tj_|_i) 


Zn  (tj+l/ti)} 


(^i+l)  {zi+1  h  [x^  (tj_|_l/tj)  , 


H  (ti+i)  dx  (f.^^) 


(2.31) 


where 

is  the  updated  state  perturbation  estimate 
is  the  propagated  state  perturbation  estimate  (dx 
explanation  following  Equation  (2.29)) 

K  (ti+i)  is  the  Kalman  filter  gain  matrix  at  time  tj+i 

Zj_|_i  is  the  measurement  realization  at  time 

h  [•,  •]  is  the  general  measurement  model  from  Equation  (2.17) 


0  from  the 


Note  that  K  (tj+i)  is  computed  from  P  and  H  (ti+i)  evaluated  along  x„  {t/U)  on  the 

interval  A  specific  equation  for  K  (ti+i)  is  given  at  the  end  of  this  subsection. 

With  the  perturbation  estimate  updated,  the  next  step  is  to  update  the  whole  state 
estimate.  As  in  the  linearized  Kalman  filter,  it  is  assumed  that  an  adequate  model  for  the 
whole  state  is  given  by 


x(t)  =  yLn{t/ti)  +  Sx.{t)  (2.32) 

Using  Equations  (2.31)  and  (2.32),  the  optimal  estimate  of  the  whole  state  is  defined  by 

x(t/U)  =  x„(t/U)  dx(t/U)  (2.33) 
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Next,  the  whole  state  estimate  is  propagated  to  the  next  measurement  time.  Since  it 
was  previously  shown  that  Sx(t/ti)  =  0  over  the  entire  interval  between  measurements, 
the  best  state  estimate  between  measurements  is  given  by  the  solution  to  Equation  (2.27) 
with  appropriate  substitutions  of  state  estimates 


x(t/ti)  =  f  [x(t/tj)  ,u(t)  ,t] 
for 

X  (ti/ti)  =  X  [tf )  (2.34) 


where 

is  the  whole  state  estimate  on  the  interval  after  measurement 

incorporation  at  time  ti 
f  [•,  •,  •]  is  the  general  dynamics  model 

X  {tf )  is  the  whole  state  estimate  immediately  following  measurement  incorporation 
at  time  U 


Finally,  the  measurement  update  at  time  tj+i  is  incorporated  using 


X 


x{ti+i/ti)  +  K  (tj+i)  {zj+i  -  h[x{ti+i/ti)  ,ti+i]} 


(2.35) 


where 

X  i4+i)  is  the  updated  whole  state  estimate  after  measurement  incorporation  at 
time  tj+i 

x(tj_|_i/ti)  is  the  propagated  whole  state  estimate  from  Equation  (2.33) 
h  [•,  •]  is  the  general  measurement  model 

The  extended  Kalman  filter  is  summarized  by  the  update  and  propagation  algorithm 
shown  in  Table  2.2. 
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Table  2.2  Extended  Kalman  Filter  Summary  [30] 


Measurement  Update 

— 

P(fjH^  [ti;S(t.  )] 

{  H  [t*;X  (t-)]  p  {t-)  (t-)]  +  R(t*) 

= 

S  (t  J  +  K  (ti)  {zi  -  h  [X  (t .  )  ,  ti]  } 

p 

= 

P(t-)  -  K(t,)H[f,;S(t-)]P(t-) 

where 

'  ' 

A 

9h[x,q] 

x=x(q-) 

State  Propagation 

Integrate  over  U  <t  <  tj+i 

x{t/ti) 

= 

{[x{t/ti)  ,u{t)  ,t] 

P  (t/ti) 

= 

F[t]x{t/ti)]P  (t/ti)  -h  P  (t/ti)  F'^  [t]x  (t/ti)] 

+G  (t)  Q  (t)  (t) 

subject  to 

^{ti/ti) 

= 

S{(+) 

P  (ti/ti) 

= 

p  K) 

where 

F[t;x(t/ti)] 

A 

5f[x,  u(^),^] 

After  Integration  Define 

s  (*i+i) 

A 

p  (*.“«) 

A 

P  {ti+l/ti) 

2.3  Multiple  Model  Adaptive  Estimator 

The  topic  of  the  Multiple  Model  Adaptive  Estimator  is  introduced  in  this  section. 
The  MMAE  was  initially  proposed  for  GPS  carrier-phase  integer  ambiguity  resolution  by 
[21]  with  a  modified  version  tested  by  [18].  The  idea  of  using  the  MMAE  (then  termed  a 
Magill  filter)  for  GPS  carrier-phase  integer  ambiguity  resolution  was  suggested  by  [21]  and 
simulations  were  performed  showing  the  technique  had  potential  for  application  in  static 
scenarios.  The  technique  was  later  modified  slightly  by  [18],  and  applied  to  kinematic  GPS 
data,  further  demonstrating  the  technique’s  potential.  This  section  provides  background 
on  the  traditional  MMAE  technique.  A  later  section  will  describe  more  fully  the  approaches 
of  [18,  21]. 
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Multiple  Model  Adaptive  Estimation  was  originally  proposed  in  [26].  The  technique 
has  since  been  described  in  a  number  of  textbooks  on  Kalman  filtering.  The  presentation 
here  is  similar  to  that  found  in  [30]. 

Consider  that  the  model  of  Equation  (2.1)  contains  an  unknown  parameter  which 
varies  slowly,  if  at  all,  when  compared  to  the  states.  Eurther,  assume  that  the  parameter  can 
take  on  values  that  are  discrete  members  of  a  finite  set  {ai,  a2, . . . ,  a^^}.  To  achieve  the  best 
state  estimation  performance,  the  unknown  parameter  must  be  estimated.  One  approach 
is  to  use  the  MMAE  shown  in  Eigure  1.2.  The  MMAE  is  a  bank  of  K  linear  Kalman  filters, 
each  built  around  a  single  member  of  the  parameter  set  {ai,a2, . . .  ,ax}.  Each  individual 
filter  is  referred  to  as  an  elemental  filter.  The  elemental  filter  outputs  are  combined  in  an 
optimally  weighted  fashion  to  produce  the  state  vector  estimate.  Additionally,  the  optimal 
weights  can  be  used  to  determine  the  best  estimate  of  the  correct  parameter  a^,  where 
k  G  {1,  2, . . . ,  K}.  With  this  introduction  in  mind,  the  MMAE  algorithm  is  derived. 

Let  a  represent  an  uncertain  parameter  in  the  system  model  of  Equation  (2.1).  In 
general,  a  may  affect  any  or  all  of  the  model  matrices:  B(^,  H,  or  R.  A  goal  of 

the  MMAE  is  to  find  the  mean  of  the  conditional  PDE  given  on  the  left  hand  side  of 

f x(ti) ,a\Z{ti)  (^)  CtjZ  (fj)  Zj)  fx(ti)\a,Z{ti)  (Cl^  ^  (j'i')  fa\Z{ti)  (^1^  (^*)  ^*) 

(2.36) 

where  /x(q),a|Z(q)  (4)C>:jZ(tj)  =  Zj)  is  the  joint  conditional  PDE  for  the  state  vector  and 
parameter  set,  x(tj)  and  a,  conditioned  on  the  realized  measurement  history  Zj. 

Equation  (2.36)  is  a  statement  of  Baye’s  rule  for  this  conditional  PDF.  In  defin¬ 
ing  the  desired  probability  density,  consider  the  first  density  on  the  right  hand  side  of 
Equation  (2.36)  /x(ti)|a,z(q)  (^l^  =  =  Zj).  Since  this  density  is  conditioned  on  a 

particular  realization  of  a,  it  is  analogous  to  the  standard  conditional  density  used  in 
Kalman  filtering  for  state  estimation.  Therefore,  under  assumptions  of  Equations  (2.1) 
through  (2.3),  it  is  a  Gaussian  density  with  mean  and  covariance  Pfc(t)'')  as  calcu¬ 

lated  by  a  standard  Kalman  filter  for  a  given  realization  of  the  parameter  a  (where  the 
subscript  k  indicates  the  dependence  upon  the  realized  value  of  the  parameter  a). 
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Next,  consider  the  second  density  on  the  right  hand  side  of  Equation  (2.36).  Although 
it  has  already  been  assumed  that  the  elements  of  a  take  on  discrete  values,  consider  for  a 
moment  the  case  in  which  a  is  a  continuous  (rather  than  discrete)  random  variable.  Using 
Baye’s  rule,  the  second  density  on  the  right  hand  side  of  Equation  (2.36)  can  be  expressed 
as 


/a|Z(q)  (Q!|Z(ti)  =  ^i)  = 

/z(q)|a,z(ti_i)  (Cl^  =  a,  Z  (ti-i)  =  Zi-i)  /a|z(q_i)  (o;|Z  (ti-i)  =  Zi-i) 

S  /z(q)|a,z(q_i)  (Cl^  Z  (tj— i)  Zi—i)  /a|z(ti_i)  (^|Z  (ii—i)  da 

A 

(2.37) 


From  the  system  model  and  assumptions  in  Equations  (2.1)  to  (2.3),  the  density 
/z(ti)|a,z(ti_i)  (Ci|3  =  CK,  Z  (tj_i)  =  Zi-i)  is  Gaussian  [29].  Its  mean  is  the  measurement 
prediction  H  (tj)  x  (t)~) ,  and  its  covariance  is  [H  (ti)  P  (tj  )  (ti)  +  R  (fi)] ,  as  computed 
in  a  Kalman  filter  based  on  the  parameter  vector  value  a.  Given  an  a  priori  density  fa  (a), 
a  state  estimate  is  given  by  the  conditional  mean  E  {x  (U)  |Z  (tj)  =  Zj}. 

Now,  consider  the  case  of  a  taking  on  discrete  values.  The  a  priori  density  for  a  can 
be  defined  by 

K 

fa(to){o^)  =  ^Pk{to)S{a  -  ak)  (2.38) 

k=l 

In  Equation  (2.38),  Pk{to)  is  the  probability  that  a  takes  on  discrete  value  at  initial 
time  to,  and  6  {a  —  a^)  is  the  Dirac  delta  function.  Thus,  fa  (a)  is  a  summation  of  Dirac 
delta  functions  weighted  by  appropriate  probabilities,  forming  a  discrete  PDF.  In  order  for 
Equation  (2.38)  to  be  a  valid  PDF,  it  must  meet  the  conditions 


Pk  {to)  >  0  y  k 

K 

E  Pk  {to)  =  1 

k=l 


(2.39) 
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Note  that  the  values  of  pk  (to)  in  Equation  (2.39)  represent  the  initial  probabilities  of  each 
element  of  the  set  {ai,a2, . . .  ,ax}-  In  the  absence  of  scientific  knowledge  or  experimen¬ 
tal  data,  the  parameter  hypotheses  are  initially  assumed  to  be  equally  likely,  leading  to 
Pk  (to)  =  M  k. 

Next,  consider  the  measurement  effects,  and  define  the  hypothesis  conditional  prob¬ 
ability 


A 


Pk{ti)  =  prob{a  =  afc|Z  (tj)  =  Zj} 


(2.40) 


Thus,  Pk  (tj)  is  the  conditional  probability  that  a  =  a^,  given  the  realized  measurement 
history  Z*. 

Substituting  Equations  (2.38)  and  (2.40)  into  Equation  (2.37),  and  using  the  sifting 
property  of  the  Dirac  delta  function,  yields 

/z(q)|a,Z(ti_i)  (^j|^  ^  {ti—l)  '^i—l)Pk 


Pk{ti)  = 


K 


(2.41) 


X]  /z(ti)|a,Z(ii_i)  (z^il^  —  3-j)Z(tj_i)  —  Zj_i)pj  (ti-l) 

i=i 


Recall  from  the  definition  of  Equation  (2.40)  that  pk  {U)  is  the  conditional  probability 
for  the  parameter  a^  given  the  measurement  history  Zj.  Thus,  if  desired.  Equation  (2.41) 
is  used  to  calculate  the  best  estimate  of  the  parameter.  The  formula  for  the  optimal 
parameter  estimate  will  be  shown  later  in  Equation  (2.48). 

Now,  consider  state  estimation.  The  optimal  state  estimate  is  derived  as  follows  [30] 


X 


(t+)  =E{x{U)\Z{U)  =  zo 

oo 

=  y /x(ti),a|Z(q)  (^,a|Z(ti)  =  Zi)dQ; 


—  OO  LA 
oo  I 


'OO 

oo 


y  ^  y  /x(ti)|a,z(q)  (4|a  =  a,z  {U)  =  Z*)  /a|z(q)  (a|Z  (U)  =  Zi)  da 


/  /x(ti)|a,z(q)  (^|a  =  a,z  (U)  =  Z*)  '^Pk  (U)  S  {a 
A  ^==1 


—  SLk)  da 


di 


2-18 


oo 

^{tf)  =  j  ^ 

—  oo 

where  Equation  (2.36)  was  applied  in  the  third  line  of  Equation  (2.42).  In  the  fourth  and 
fifth  lines  of  Equation  (2.42),  /a|z(ti)  (q;|Z  {U)  =  Zi)  is  replaced  by  J2k  =  i  Pk  {U)  6  {a  —  a^) 
by  invoking  the  definition  of  Pk{ti)-  Also,  rather  than  writing  a  series  of  integrals  over 
each  element  of  the  vector  the  common  engineering  notation  of  writing  a  single  integral 
with  boldface  limits  (  f  {•}  )  was  used.  Assuming  ^  is  n-dimensional,  this  notation  is 

—  OO  ^ 

interpreted  as 


K 

E 

.k=l 


/x(ti)|a,Z(ti)  (^1^  ^kj  Z  (tj)  Zi)  Pk  {ti) 


di 


(2.42) 


oo  oo  oo  oo 

I  {■}d^  =  I  I  ...  I  {■}dCldC2.-.dCn 

—  OO  — oo  — oo  — oo 

Since  the  summation  in  Equation  (2.42)  is  finite,  the  order  of  integration  and  sum¬ 
mation  can  be  interchanged  [1]  yielding 


K 


s(4)  =  E 

k=l 


f  ^/x(q)|a,z(q)  (C|a  =  ^k,  Z  (ti)  =  Zi) 


Pk  (ti 


(2.43) 


The  term  inside  the  square  bracket  of  Equation  (2.43)  is  simply  the  optimal  state  estimate 
conditioned  upon  a  given  parameter  value  and  the  realized  measurement  history  Zj. 
Define  this  quantity  to  be  (tf ) .  The  updated  state  estimate  then  takes  on  its  final  form 
given  by 

K 

^{4)  =  '^^k{tf)pk(ti)  (2.44) 

k=l 

Equation  (2.44)  shows  that  the  final  updated  state  estimate  is  the  weighted  sum  of 
the  elemental  filter  state  estimates,  weighted  by  the  probability  that  the  parameter  for 
that  particular  elemental  filter  is  the  correct  parameter  (based  on  the  measurements  seen 
to  date).  This  is  exactly  the  form  that  is  shown  in  Eigure  1.2. 
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Next,  consider  the  covariance  update  for  the  MMAE.  From  the  definition,  the  co¬ 


variance  is 


P  (4)  =  -x(t+)]  [x(ti)  -S(t+)]^|Z(ti)  =  Zij 

oo 

=  /  [£-S(^)]  [S-SW)]^/x«.)|Z(y(«|Z(t.)  =  Zi)rf« 

—  OO 

oo 

=  /  [^-x(t+)]  [^-X(t+)]^y/x(t,),a|Z{ti)(^,«|Z(ti)  =Zi)dQ 

-oo  A 

oo 

=  /  [«  -  S  (<.*)]  [f  -  s  (*.")]’■ 

—  oo 

fx(ti)\a,Z{ti)  (^1^  Z  (tj)  Zj)  /a|Z(ti)  (*^|Z  (tj)  Zj)  dOL  d^ 
A 


K 


=  I  ^Liu)\a,ziu){^\a  =  ak,Z{ti)  =  Zi)Pk{ti)  d$, 

(2.45) 


k=l 


Again,  since  the  summation  is  finite,  the  order  of  integration  and  summation  can  be 
interchanged.  Additionally,  /a|z(ti)  (c>:|Z(tj)  =  Zj)  may  be  brought  outside  the  integral  in 
the  form  pk  (U)  using  Equation  (2.40).  Equation  (2.45)  then  becomes 

^  (  7  1 

^(7)  =  '^Pk{ti)l  /  f^(^ti)\a,z{u){^\a  =  a.k,Z{ti)  =  Zi)d^' 

fc=l  Uoo  J 

K  f  °° 

=  ^Pk  {ti)  I  j  ^^^/x(ti)|a,z(q)  (^|a  =  afc,  Z  {ti)  =  Zi)  d$ 

V— OO 
oo 

~  J  ^  )  ^^/x(q)|a,z(q)  (C|a  =  afc,  Z  (U)  =  Zi)  d$ 

—  OO 
OO 

-  j  {7 )  /x(ti)|a,z(q)  (^|a  =  afc,  Z  {U)  =  Zi)  d$ 

—  OO 

+  y  X  {7 )  {4 )  /x(q)|a,Z(q)  (^|a  =  afc,  Z  (t*)  =  Zi)  d$  | 
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K 


^(4)  =^Pk{ti)l  I  C^'^/x(ti)|a,Z(ti)(^|a  =afc,Z(ti)  =  Zi)d^ 

V— oo 
oo 

“  ^  {4 )  j  fx{ti)\a,z(ti)  (^|a  =  afc,  Z  (ti)  =  Zi) 

—  OO 

oo 

J  ^Mti)\a,ziti)  (C|a  =  afc,  Z  {ti)  =  Zi)  Sf'  {tf ) 

-OO 

OO 

+  X  {tf )  {tf )  J  f„(ti)\a,z{ti)  (C|a  =  afc,  Z  (ti)  =  Zi)  d^ 

—  OO  ^ 

K  (  ^ 

=  '^Pk  {ti)  I  j  ^^'^/x(ti)|a,z(ti)  (^|a  =  afc,  z  (ti)  =  Zi)  -  X  (t+)  (tf) 

V~oo 

-Sfc  (t+ )  (t+ )  +  X  {tf )  X^  (t+ )  I 

P  )  =  '^Pk  (**)  {Pfc  )  +  [xfc  )  -  S  [Sfc  )  -  X  )] ^}  (2.46) 


fc=i 


Thus,  the  state  estimate  error  covariance  of  Equation  (2.46)  is  dependent  upon  the 
overall  state  estimate  X  (^tf ) ,  the  state  estimate  from  each  elemental  filter  X^  (^tf ) ,  the 
error  covariance  from  each  elemental  filter  {tf  ),  and  the  probability  pk  {tf  that  each 
elemental  filter  contains  the  correct  parameter  value.  The  estimates  X^  {tf)  and  {tf) 
are  readily  available  from  each  elemental  Kalman  filter.  If  pk  {tf  is  available,  then  the  filter 
shown  in  Figure  1.2  is  defined.  Note  that  P  {tf)  is  not  required  for  an  online  MMAE.  As 
a  final  step,  consider  the  calculation  of  pk  {tf. 

Given  an  initial  estimate  pk  (to)  subject  to  conditions  of  Equation  (2.39),  the  proba¬ 
bilities  Pk  {tf  can  be  updated  at  each  measurement  epoch  using  Equation  (2.41).  The  only 
value  still  needed  for  the  update  is  the  density  /z(q)|a,z(q_i)  (zi|a  =  a^,  Z  {U-f  =  Zj_i)  for 
k  =  1,  2,  ...  K,  which  is  calculated  as  [30] 
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(2.47) 


fz(ti)\a,Z(ti-i)  Zj_i) 

=  |i  exp  I  -  (t- )  (t- )  rfc  (t” )  | 

(27r)2  |Afc(tJ  |2  I  2  J 

where  A^  (t“ )  =  (tj)  (t“ )  (ti)  +  R  (ti)  is  the  covariance  matrix  of  the  mea¬ 
surement  residual  vector  (t“)  from  elemental  filter  k.  Both  (t“)  and  A^  {t~)  are 
calculated  from  each  of  the  K  elemental  filters,  and  (t“)  =  Zj  —  H  (tj)  (t~),  where 
Zj  is  the  realized  measurement  vector  at  ti. 

If  desired,  the  best  parameter  estimate  can  be  derived  as  follows: 


A 


a{ti)  =  E{a{ti)\Z{U)  =  Z,} 

oo 

f  *^/a|Z(ti) 


-oo 

oo  r  K 


—  oo 
K 


'^Pk{ti)6{a  - 
.fc=i 

a(ti)  =  '^akPk{ti) 


da 


k=l 


(2.48) 


As  in  previous  derivations,  the  last  line  of  Equation  (2.48)  is  the  result  of  interchanging 
the  order  of  integration  and  summation,  and  applying  the  sifting  property  of  the  Dirac 
delta  function. 

Finally,  the  parameter  estimate  precision  is  given  by  the  parameter  conditional  co- 
variance 


Pa(ti)  =  E|[a  -a(ti)]  [a  -  a  (t*)]'^  |Z  (tj)  =  Z^j 

=  E  [aa^|Z  (ti)  =  Z*]  -  a  {U)  E  [a^|Z  {U)  =  Z*] 
-  E  [a|Z  (ti)  =  Zi]  a'^  {U)  +  a  (U)  a'^  (U) 
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K  K 

Pa  {ti)  =  ^  akalpk  (ti)  -  a  {U)  ^  alpk  (t*) 
k=l  k=l 

K  K 

-  ^  akPk  (ti)  (ti)  +  a  (U)  a^  (U)  '^Pk  (U) 
k=l  k=l 

K 

=  ^  [a^a^  -  a  {U)  al  -  a^a^  {U)  +  a  {U)  (t*)]  pk  {U) 

k=l 
K 

Pa(t*)  =  '^[ak-a^  {ti)][ak-a^  {ti)]^  Pk{ti)  (2.49) 

k=l 

Note  that  Equations  (2.48)  and  (2.49)  need  not  be  calculated  to  estimate  the  state 
X  (t+).  Also,  care  should  be  taken  when  tuning  the  MMAE  [30].  Dynamics  pseudonoise, 
added  to  account  for  modeling  inadequacies,  has  a  tendency  to  make  measurement  residuals 
from  elemental  filters  have  similar  magnitudes.  However,  the  MMAE  relies  on  the  weighted 
sum  of  squared  residuals  (also  called  the  likelihood  quotient)  (t“)  being 

substantially  different  between  the  elemental  filter  based  on  the  correct  hypothesis  and 
any  other  elemental  filter.  Excessive  amounts  of  dynamics  pseudonoise  being  added  can 
incapacitate  the  adaptivity  of  the  MMAE.  So,  “conservative”  tuning  (adding  sufficient 
pseudonoise  to  the  model  to  preclude  filter  divergence)  is  to  be  avoided.  Instead,  “tight” 
filter  tuning  without  much  (if  any)  added  pseudonoise  is  preferred.  Therefore,  to  prevent 
the  MMAE  from  overemphasizing  incorrect  elemental  filters,  care  must  be  taken  when 
tuning  the  MMAE. 

Einally,  in  this  derivation,  it  was  assumed  that  the  correct  parameter  is  a  member 
of  a  discrete  set  of  parameters  (i.e.,  a  G  {ai,a2, . . .  ,ai^}).  When  this  assumption  is  met, 
the  MMAE  will  converge  to  the  correct  parameter  value  a^.  However,  if  this  assumption 
is  not  met,  the  MMAE  will  converge  to  the  parameter  value  which  is  closest  to  the  correct 
value  in  the  Baram  distance  measure  sense  [5,  6]. 


2-4  MMAE  Performance  Enhancements 

Several  modifications  to  the  traditional  MMAE  have  been  proposed  to  enhance  per¬ 
formance.  Most  are  ad  hoc  methods  which  are  devised  to  address  problems  in  specific 
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applications.  Good  summaries  of  these  enhancements  are  found  in  [31,  46].  Some  of  these 
enhancements  are  listed  and  discussed  next. 

1.  f3  Stripping 

2.  Probability  Lower  Bounding 

3.  Kalman  Filter  Tuning 

4.  Elemental  Filter  Gain  Modulation 

5.  Scalar  Penalty  Increase 

6.  Probability  Smoothing 

7.  Increased  Residual  Propagation 

8.  Interacting  Multiple  Model 

2.4-1  /3  Stripping.  While  using  the  MMAE  for  sensor  failure  identification, 

[32,  45]  found  that  the  MMAE  was  susceptible  to  false  alarms  when  the  likelihood  quotients 
from  all  the  elemental  filters  had  similar  magnitudes  over  time.  Under  this  condition, 
the  elemental  filters  are  all  equally  adequate,  and  the  MMAE  should  equally  distribute 
the  hypothesis  conditional  probabilities  among  the  filters.  However,  the  MMAE  instead 
assigns  more  probability  weight  to  the  elemental  filter  with  the  smallest  jA^  j.  This 
is  the  so-called  “/^-dominance”  effect  [2,  3,  32]  caused  when 

"  (27r)W2|Afc  (t-)lV2 

dominates  over 

exp{-^r^  {t~ )  A-i  {t- )  Tk  (t“)} 

in  the  hypothesis  conditional  probability  calculation.  Such  behavior  by  the  MMAE  is 
inappropriate  and  undesirable  because  the  filter-computed  residual  covariance  A^  (t“)  is 
completely  divorced  from  the  actual  sensor  measurements.  In  fact,  if  linear  Kalman  filters 
are  used  for  the  elemental  filters,  then  A^  is  pre-computable  before  any  measurements 
are  taken. 


2-24 


To  prevent  the  /3-dominance  effect,  [32,  45]  stripped  the  (3k  (U)  terms  from  all  the 
Gaussian  densities  in  Equation  (2.41).  Thus, 

exp  I  -^rl  {t- )  {t- )  Vk  (t“ )  [ 

(2vr)  2|Afc  (t,  )  h  2  J 

was  replaced  by 

exp  |-^r^  {t- )  A-i  {t- )  Vk  (t“ )  | 

in  the  Gaussian  densities  of  Equation  (2.41).  While  these  stripped  “densities”  do  not 
integrate  to  1,  and  therefore  are  not  true  probability  densities,  the  modified  hypothesis 
conditional  probability  formula  is  still  valid  because  the  denominator  still  scales  the  prob¬ 
abilities  associated  with  each  elemental  filter  such  that  they  all  sum  to  1. 

As  mentioned  at  the  beginning  of  this  section,  sensor  failure  identification  is  an  appli¬ 
cation  in  which  /3-dominance  has  been  noted  and  /3-stripping  applied.  In  that  application, 
failed  sensors  are  commonly  modelled  in  the  elemental  filters  by  setting  all  elements  in 
the  appropriate  row  of  the  measurement  matrix  Hj  [ti)  to  zero.  This  has  a  tendency  to 
make  the  determinant  of  the  filter-computed  residual  covariance  matrix  associated  with 
that  filter,  |Aj  |,  smaller  than  the  same  quantity  for  the  fully  functional  measurement 
model.  Thus,  for  this  application,  the  MMAE  tends  to  be  pre-disposed  to  declare  sensor 
failures  even  if  none  exist,  due  solely  to  /3-dominance.  This  effect  was  noted  in  an  MMAE 
used  for  flight  control  sensor  failures  [32]. 

Although  (3  stripping  is  not  used  in  this  research,  the  concept  of  /3-dominance  is 
discussed  frequently.  This  research  shows  that  the  7  design  parameter  in  the  GRMMAE 
derived  in  Ghapter  III  causes  the  GRMMAE  to  be  highly  susceptible  to  /3-dominance.  This 
fact  is  derived,  explained,  and  illustrated  in  Ghapters  III  and  IV. 

2.4-2  Probability  Lower  Bounding.  Observe  in  Equation  (2.41)  that,  if  PkiU) 
goes  to  zero,  it  will  remain  zero  for  all  t  >  ti.  This  causes  elemental  filter  k  to  be  removed 
from  the  MMAE.  To  prevent  this  from  happening,  it  is  common  practice  to  set  a  lower 
bound  e  on  pk  (U)  (typically  0.001  <  e  <  0.01).  If  a  calculated  value  of  pk  (U)  is  less  than 
e,  then  pk  (ti)  is  set  equal  to  e,  and  all  remaining  probabilities  are  rescaled  to  maintain 
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J2k=iPk  i^i)  ~  Probability  lower  bounding  is  implemented  only  in  Section  4.2.  Aside 
from  that  example,  probability  lower  bounding  is  not  used  in  this  research. 

2.4-3  Kalman  Filter  Tuning.  As  mentioned  in  Section  2.3,  the  elemental  filters 
must  be  carefully  tuned.  Addition  of  too  much  process  pseudonoise  to  account  for  modeling 
inadequacies  can  have  the  effect  of  making  the  likelihood  quotients  from  all  elemental  filters 
appear  similar  over  time.  This  degrades  the  MMAE’s  ability  to  distinguish  between  models 
and  converge  to  the  correct  hypothesis.  In  the  failure  detection  and  isolation  application, 
this  problem  manifests  itself  as  detection  delays  or  missed  detections.  However,  properly 
tuning  the  elemental  filters  reduces  false  and  missed  alarms  [8,  9,  16,  27,  30,  28].  Generally, 
designers  use  ad  hoc  methods  to  determine  appropriate  elemental  filter  noise  covariance 
matrices  based  on  physical  insight  into  the  system.  Assuming  elemental  filter  k  is  correct, 
the  goal  in  tuning  the  elemental  filters  is  to  choose  the  noise  covariance  matrices  such 
that  rl  {t- )  A-i  {t- )  Yk  {t- )  is  on  the  order  of  the  measurement  dimension  m,  while  the 
likelihood  quotients  from  the  other  elemental  filters  are  larger  than  m  [46].  Kalman  filter 
tuning  is  used  in  this  research. 

2.4.4  Elemental  Filter  Gain  Modulation.  Related  to  the  previous  issue  of  elemen¬ 
tal  filter  tuning,  [25]  derived  a  method  for  modulating  the  Kalman  gain  in  the  elemental 
filters  to  increase  hypothesis  distinguishability.  Assume  that  the  parameter  of  interest 
affects  the  elemental  filter  dynamics  model.  Then,  to  increase  the  MMAE’s  ability  to  dis¬ 
tinguish  differences  between  the  elemental  filters,  one  wants  to  adjust  the  Kalman  gain  in 
the  elemental  filters  to  emphasize  the  dynamics  model  and  de-emphasize  the  measurement 
update.  This  is  accomplished  by  setting  the  dynamics  driving  noise  covariance  matrix 
Qd  {ti)  artificially  low.  However,  this  has  the  drawback  of  overemphasizing  the  dynamics 
model  during  the  state  estimate  measurement  update,  degrading  state  estimation  perfor¬ 
mance. 

In  this  method,  the  Kullback  Discrimination  Index  (KDI)  (similar  to  a  distance 
metric,  but  not  a  true  metric)  is  used  after  each  measurement  update  to  measure  the 
“distance”  between  elemental  filters  based  on  the  likelihood  quotients  from  the  elemental 
filters.  Assuming  that  adequate  state  estimation  performance  is  being  achieved,  if  the 
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elemental  filters  are  deemed  “too  close”  for  good  hypothesis  distinguishability  (i.e.,  the 
likelihood  quotients  of  two  or  more  elemental  filters  are  of  similar  magnitudes),  then  the 
dynamics  noise  covariance  matrices  in  the  elemental  filters  are  decreased  for  the  next 
measurement  update.  This  has  the  effect  of  increasing  the  emphasis  on  the  elemental  filter 
dynamics  model  during  the  measurement  update.  This  increased  model  emphasis  increases 
the  MMAE’s  ability  to  distinguish  between  elemental  models  (hypotheses).  At  the  next 
measurement  update,  the  KDI  is  calculated  again  and  the  process  is  repeated.  During  this 
process  of  decreasing  the  dynamics  pseudonoise,  care  must  be  taken  to  ensure  that  Qd  (U) 
is  not  decreased  to  the  point  where  adequate  state  estimation  performance  is  no  longer 
achieved.  Therefore,  this  method  modulates  the  elemental  filter  Kalman  gain  at  each 
measurement  to  maximize  hypothesis  distinguishability  while  ensuring  state  estimation 
performance  is  not  degraded.  Elemental  filter  gain  modulation  is  not  used  in  this  research. 

2.4-5  Scalar  Penalty  Increase.  A  method  for  speeding  MMAE  convergence,  sug¬ 
gested  by  [31],  is  scalar  penalty  increase.  In  this  method,  exp{— {t~)} 
in  the  Gaussian  densities  of  Equation  (2.41)  is  replaced  by  exp{  —ar^  (t~)  (t~)  {t~)} 

where  a  is  some  scalar  larger  than  1/2.  Similar  to  the  /3-stripped  “densities”  in  Section 
2.4.1,  these  “densities”  with  increased  scalar  penalties  are  not  proper  probability  densities 
(i.e.  do  not  integrate  to  1).  However,  due  to  the  normalization  terms  in  the  denominator 
of  Equation  (2.41),  the  pk  {ii)  values  will  still  sum  to  1  in  the  MMAE. 

Making  a  greater  than  1/2  assigns  higher  penalties  to  large  likelihood  quotients, 
compared  with  a  =  1/2.  While  this  speeds  convergence  to  the  correct  hypothesis,  it  can 
also  cause  large  conditional  probability  fluctuations  during  convergence.  If  the  MMAE  is 
used  for  failure  detection,  these  fluctuations  may  be  interpreted  as  momentary  false  alarms. 
Thus,  the  a  value  (which  is  chosen  through  ad  hoc  methods)  should  be  chosen  carefully 
to  meet  some  desired  performance  specification.  Scalar  penalty  increase  is  an  alternative 
to  GRMMAE  derived  in  this  research.  It  can  achieve  faster  convergence  without  the 
susceptibility  to  /3-dominance  which  plagues  the  GRMMAE.  However,  care  must  be  taken 
to  avoid  the  false  alarm  problem  induced  by  scalar  penalty  increase.  An  example  of  scalar 
penalty  increase  is  given  in  Section  4.5. 
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2.4-6  Probability  Smoothing.  This  method  is  used  primarily  for  failure  detection. 
After  a  system  changes  operating  modes,  the  hypothesis  conditional  probabilities  for  the 
elemental  filters  will  undergo  transients.  These  transients  may  cause  false  alarms  [15,  31, 
33,  34,  43,  44].  To  decrease  these  false  alarms,  the  hypothesis  conditional  probabilities  can 
be  smoothed  using  a  moving  window  average.  The  window  size  is  generally  determined 
empirically.  As  the  window  size  increases  (i.e.,  more  samples  are  used  to  calculate  the 
average),  the  resistance  to  false  alarms  increases.  However,  this  also  causes  slower  response 
to  actual  failures.  This  research  does  not  use  probability  smoothing. 

2.4.7  Increased  Residual  Propagation.  Increased  residual  propagation  [15,  31] 
is  similar  to,  yet  less  complex  than,  elemental  filter  gain  modulation.  In  this  method, 
the  elemental  filter  state  estimates  are  propagated  a  few  sample  times  without  allowing 
measurement  updates.  By  doing  this,  the  likelihood  quotients  from  incorrect  elemental 
filters  grow  because  the  measurement  updates  are  not  correcting  the  state  estimates  toward 
the  measurements.  However,  the  likelihood  quotient  from  the  correct  elemental  filter  will 
remain  relatively  small,  because  its  hypothesized  model  matches  the  measurements.  Thus, 
hypothesis  distinguishability  is  increased,  leading  to  faster  convergence.  The  number  of 
samples  to  propagate  before  measurement  incorporation  is  determined  empirically.  The 
drawback  to  this  method  is  that  it  increases  hypothesis  conditional  probability  fluctuations, 
which  may  lead  to  false  alarms  in  failure  detection  applications  [15,  31].  Additionally,  state 
estimation  performance  will  decrease  due  to  the  increased  time  between  measurement 
updates.  This  research  does  not  use  increased  residual  propagation. 

2.4.S  Interacting  Multiple  Model.  The  MMAE  assumes  that  the  unknown  pa¬ 
rameter  a  varies  slowly,  if  at  all,  compared  to  the  states.  A  similar  estimator  exists  which 
allows  the  parameter  to  vary  under  certain  assumptions.  This  estimator  is  the  Interacting 
Multiple  Model  (IMM)  [4]. 

The  IMM  is  similar  to  the  MMAE  in  the  sense  that  it  uses  multiple  parallel  Kalman 
filters  each  built  around  a  potential  parameter  realization  a^.  Unlike  the  MMAE,  however, 
the  IMM  assumes  that  the  hypothesis  conditional  probability  vector  (composed  of  the 
probabilities  associated  with  each  elemental  filter)  propagates  according  to  a  Markovian 
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model.  The  IMM  recursively  produces  updated  state  and  covariance  estimates  based  on  the 
previous  state  and  covariance  estimates  and  the  previous  parameter  estimate.  However, 
IMM  is  suboptimal,  because  it  only  produces  estimates  conditioned  on  the  current  and 
most  recent  past  moment  in  time,  rather  than  the  entire  past  (see  [4]  for  an  explanation). 

This  dissertation  was  inspired  by  the  GPS  carrier-phase  integer  ambiguity  resolution 
application.  In  that  application,  the  parameter  is  the  GPS  carrier-phase  integer  ambiguity 
vector  which  is  constant.  While  the  IMM  allows  added  flexibility  compared  to  the  MMAE, 
it  is  inappropriate  for  that  application  because  the  parameter  is  truly  constant.  Although 
the  GRMMAE  is  general  enough  to  be  applied  to  other  applications,  because  it  was  inspired 
by  GPS  carrier-phase  integer  ambiguity  resolution,  this  research  does  not  apply  the  IMM. 

2.5  Global  Positioning  System 

GPS  is  a  direct  sequence  spread  spectrum  (DSSS)  satellite  navigation  system  con¬ 
sisting  of  24  or  more  medium  earth  orbiting  satellites  providing  navigation  information  to 
passive  receivers.  The  system  operates  using  the  principle  of  trilateration.  The  satellite 
embeds  precise  timing  information  within  the  transmitted  signal.  The  receiver  uses  signal 
propagation  delay  to  calculate  the  satellite  distance  using  the  formula  r  =  Vp{tt  —  U), 
where  r  is  the  range  from  the  receiver  to  the  satellite,  Vp  is  the  velocity  of  propagation 
(usually  assumed  to  be  the  speed  of  light  in  a  vacuum),  tt  is  the  transmit  time,  and  tr  is 
the  reception  time.  In  addition  to  transmitting  precise  timing  information,  each  satellite 
transmits  ephemeris  information  which  is  used  to  calculate  the  satellite’s  position.  Prom 
knowledge  of  the  positions  of  at  least  four  satellites  and  the  range  to  each  satellite,  a 
three-dimensional  user  position  plus  clock  bias  solution  is  estimated. 

2.5.1  The  Global  Positioning  System  Signal.  Before  examining  the  GPS  signal, 
consider  a  general  DSSS  signal,  represented  by 

s{t)  =  \/2Pcos\LVQt  +  9c{t)  -|-  9(i{t)\  (2.50) 

where  s{t)  is  the  transmitted  signal,  P  is  the  average  carrier  power,  lvq  is  the  carrier 
angular  frequency,  9c  (t)  is  the  spreading  code  waveform,  and  9d  (t)  is  the  data  modulation. 
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Table  2.3  GPS  Spreading  Code  Characteristics  [22] 


Parameter 

C/A  Code 

P  Code 

Data  Rate 

50  Hz 

50  Hz 

Chip  Rate 

1.023  Mchips/sec 

10.23  Mchips/sec 

Code  Period 

1023  Chips  (1  ms) 

6  X  10^^  Chips  (1  week) 

Code  Chip  Length 

293.25  m 

29.325  m 

Carrier  Band  Designation 
and  Wavelength 

LI 

1575.42  MHz 
0.1904  m 

LI,  L2 

1575.42,  1227.6  MHz 
0.1904,  0.2444  m 
respectively 

Assuming  antipodal  binary  phase  shift  keying  (BPSK)  is  used  for  the  spreading  and  data 
modulations,  an  equivalent  form  of  Equation  (2.50)  is  given  by 


s{t)  =  \/2Pc{t)d{t)cos{ujQt) 
c{t)  =  ±  1 

d{t)  =  ±1  (2.51) 

where  d{t)  and  c{t)  are  the  data  and  spreading  waveforms,  respectively. 

The  GPS  consists  of  satellites  transmitting  two  separate  DSSS  signals  of  the  form 
of  Equation  (2.51).  The  first  signal,  called  the  Precision  code  (P-code),  is  generally  en¬ 
crypted^  and  reserved  for  military  users.  The  second,  the  Coarse/ Acquistion  (C/A)  code, 
is  available  to  civilian  users.  The  P-code  is  transmitted  at  two  carrier  frequencies,  LI  and 
L2,  and  C/A  code  is  transmitted  at  LI.  The  GPS  signal  properties  are  summarized  in 
Table  2.3. 

2.6  Precise  GPS  Positioning 

GPS  positions  can  be  calculated  using  either  the  code  phase  or  the  carrier-phase  for 
timing.  The  code  phase  (or  pseudorange^)  measurement  is  used  for  general  navigation, 

^Encrypted  P-code  is  called  “Y-code” 

GPS  range  measurement  is  termed  a  “psendorange”  measurement  because  it  is  the  range  to  the 
satellite  incorporating  the  effect  of  a  clock  bias. 
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Table  2.4  GPS  Code  Phase  and  Carrier-Phase  Positioning  Accuracy[18] 


Standard 
Code  Phase 
GPS 

Code  Phase 
Differential  GPS 

Carrier-Phase 
Differential  GPS 
Fixed  Ambiguities 

Approximate 
3-D  Accuracy 

6-10  m 

1  m  [36] 

0.01-0.03  m 

producing  position  estimates  with  errors  on  the  order  of  meters.  However,  due  to  the  rel¬ 
atively  short  carrier  wavelength,  carrier-phase  can  be  used  to  produce  position  estimates 
with  errors  on  the  order  of  tens  of  centimeters.  Note,  though,  that  carrier-phase  based 
(or  any  other  differential  GPS  based)  positions  are  positions  of  one  point  relative  to  an¬ 
other.  Table  2.4  summarizes  typical  3-D  accuracies  achieved  using  stand-alone  code  phase, 
differential  code  phase,  and  differential  carrier-phase  positioning  methods^.  Differential 
(code  or  carrier- phase)  GPS  is  a  process  of  using  multiple  receivers  and/or  satellites  to 
remove  common  errors.  Single  and  double  difference  carrier-phase  differential  GPS  will  be 
described  shortly. 

2.1  Integer  Ambiguities 

To  gain  an  understanding  of  carrier-phase  integer  ambiguities,  consider  Figure  2.1 
which  shows  a  GPS  signal  propagating  from  satellite  k  to  receiver  A.  The  range  to  the  GPS 
satellite  is  measured  by  counting  the  number  of  carrier-phase  cycles  between  the  satellite 
and  the  receiver,  and  multiplying  by  the  carrier  wavelength.  It  is  assumed  that  the  receiver 
begins  tracking  the  satellite  signal  at  time  to  and  continues  to  track  until  time  t.  Once 
tracking  is  established  at  to,  the  receiver  counts  the  carrier  cycles  received  from  that  point 
forward  in  time.  Additionally,  the  receiver  is  able  to  determine  where  in  the  partial  cycle 
it  begins  to  track,  depicted  as  (j)frac{to)  in  Figure  2.1.  However,  the  receiver  is  unable  to 
determine  the  integer  number  of  cycles  N  which  have  already  propagated  at  time  to-  This 
integer  number  of  cycles  is  the  earrier-phase  integer  ambiguity  which  must  be  resolved  to 

^Table  2.4  assumes  that  the  intentional  degradation  of  Selective  Availability  (SA)  is  not  present.  SA 
was  turned  off  in  2000  and  is  not  expected  to  be  used  again.  However,  the  capability  to  use  SA  again  does 
exist,  in  which  case  standard  C/A  code  phase  position  accuracies  would  be  on  the  order  of  100  m,  while 
the  differential  (code  and  carrier-phase)  accuracies  would  remain  the  same. 
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Figure  2.1  Carrier-Phase  Measurement  Example 


determine  range  to  the  satellite.  Once  N  is  determined,  it  will  remain  constant  as  long  as 
continuous  carrier  tracking  is  maintained. 

The  carrier-phase  is  calculated  in  receiver  A  by  integrating  the  Doppler  shift  present 
in  the  received  signal  as  described  by  Equation  (2.52)  [38], 

t 

(p{t)  =  J  Afmeas{r)dT  +  (J)  {to)  N  (2.52) 

^0 

^ - V - ^ 

4^meas  it) 

where  cj)  {t)  is  the  total  carrier- phase  cycle  count  at  time  t  (corresponding  to  the  range 
between  the  satellite  and  the  user) ,  Afmeas  (t)  is  the  received  signal  Doppler  shift  measured 
by  the  receiver  at  time  r,  (j)  (to)  is  the  carrier  cycle  count  at  to,  which  may  include  both  an 
integer  portion,  (pint  {to),  and  a  partial  cycle  portion,  (pfrac{to),  and  N  is  the  ambiguous 
integer  carrier  cycle  count  at  the  initial  time.  The  quantity  (pmeas  {t)  =  Afmeas  {t)  dr  + 
(p  {to)  is  the  carrier-phase  measurement  available  from  receiver  A.  One  application  for  the 
GRMMAE  is  resolving  the  integer  ambiguity  N  in  Equation  (2.52)  and  Eigure  2.1. 
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2.8  Differential  Carrier- Phase  GPS 

For  ambiguity  resolution  to  work  efficiently,  GPS  range  errors  must  be  removed 
or  reduced  prior  to  determining  the  carrier-phase  integer  ambiguities.  This  is  normally 
accomplished  through  differential  GPS  procedures  known  as  single  and  double  differenc¬ 
ing.  Single  differencing  is  described  first,  followed  by  a  description  of  double  differencing. 
The  concepts  are  presented  for  carrier-phase  measurements;  however,  these  techniques  are 
equally  valid,  and  similarly  derived,  for  code  phase  measurements. 

2.8.1  Single  Differenee  Carrier-Phase  Differential  GPS.  Gonsider  two  GPS  re¬ 
ceivers  tracking  a  signal  from  a  common  GPS  satellite  as  shown  in  Figure  2.2.  A  single 
difference  is  formed  by  subtracting  the  carrier-phase  measurement  of  one  receiver  from  the 
other.  Using  notation  of  [35],  Equation  (2.53)  shows  the  carrier-phase  measurement  from 
receiver  u  to  satellite  k,  in  units  of  carrier  cycles: 


Ak) 

(pu 


=  A 


-1 


(fc) 

rii  - 


f) 


(k) 


+  +  m 


(k) 

(f),U 


+  f  {6tu  -  St^^^)  + 


(2.53) 


where 

k  as  a  superscript,  indicates  the  k^^  satellite  being  tracked 

u  as  a  subscript,  indicates  the  receiver  tracking  the  satellite  {A  or  B) 

is  the  carrier-phase  multipath  tracking  error  for  receiver  u  (m) 

A  is  the  carrier  wavelength  (0.19  m  for  Li  or  0.24  m  for  L2) 

r  is  the  true  range  from  the  receiver  to  the  satellite  (m) 

I  is  the  ionospheric  delay  parameter,  empirically  determined  to  be 

40.3  TEG  (Hz2.  m)[37] 

T  is  the  tropospheric  error  in  the  measurement  (m) 

/  is  the  carrier  frequency  (Hz) 

6tu  is  the  receiver  clock  error  (s) 
is  the  satellite  clock  error  (s) 


2-33 


/ 


/  \ 
/  \ 


/ 


k  / 

/ 


/ 


/ 


/ 


\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

\ 

_ \ 


A  B 

Figure  2.2  Single  Difference  GPS  Geometry 


N  is  the  carrier-phase  integer  ambiguity 

is  the  remaining  carrier-phase  measurement  error  (measurement  noise  plus 
possible  tracking  loop  errors  under  high  dynamics)  (cycles) 

If  receivers  A  and  B  are  separated  by  a  short  baseline  (generally  on  the  order  of 
tens  of  kilometers),  many  of  the  errors  in  Equation  (2.53)  are  highly  correlated.  Thus,  by 
subtracting  the  phase  measurement  of  one  receiver  from  that  of  another,  correlated  errors 
can  be  essentially  removed.  Note  that  this  difference  measurement  provides  an  indication 
of  the  relative  position  of  the  two  receivers,  not  the  range  to  the  satellite.  Thus,  when 
differencing  is  used,  the  position  of  one  receiver  (termed  a  “rover”)  is  determined  relative 
to  another  receiver,  which  is  typically  at  a  known  location  (termed  a  “reference”).  Letting 
A  indicate  a  single  difference  between  receivers.  Equation  (2.54)  shows  the  single  difference 
measurement  between  receivers  A  and  B: 
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In  Equation  (2.54),  A(/>^^  is  the  single  difference  carrier-phase  measurement  between 
(k) 

receivers  A  and  B,  is  the  single  difference  phase  integer  ambiguity  between  receivers 

A  and  B,  Av^^{^  =  and  the  A  preceding  the  error  terms  indicates  residual 

errors  after  differencing  (i.e.,  the  single  differenced  errors).  Theoretically,  for  two  receivers 
at  the  same  location,  these  residual  errors  are  zero®.  As  the  baseline  distance  between  the 
receivers  increases,  the  error  correlation  between  the  two  receivers  decreases  and  residual 
errors  grow.  Finally,  since  AN^^  is  the  result  of  differencing  two  integers,  it  remains  an 
integer  quantity. 

2.8.2  Double  Difference  Carrier-Phase  Differential  GPS.  Equation  (2.54)  con¬ 
tains  the  receiver  clock  error  5tA  —  hts,  which  cannot  be  removed  by  single  differencing 
between  receivers.  To  obtain  the  most  accurate  relative  position,  this  error  is  normally 
removed  by  forming  a  double  difference,  the  difference  between  single  difference  measure¬ 
ments  using  two  satellites  as  shown  in  Figure  2.3.  Letting  VA  represent  a  double  difference. 
Equation  (2.55)  gives  a  mathematical  description  of  double  differencing  between  receivers 
A  and  B  and  satellites  k  and  j. 


VAf 


(kj) 

AB 


A2b  -  A2b 


®For  this  to  be  completely  correct,  the  two  receivers  would  have  to  use  the  same  antenna,  so  that  the 
multipath  error  would  be  identical  in  the  two  receivers. 
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Figure  2.3  Double  Difference  GPS  Geometry 
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(2.55) 


In  Equation  (2.55),  VA$^^^  is  the  double  difference  phase  measurement  between 

receivers  A  and  B  and  satellites  k  and  j.  Note  that,  as  in  the  single  difference  case, 

(kj) 

the  double  difference  ambiguities,  V AN^^\  are  again  integers.  Also,  the  VA  preceding 
the  error  terms  on  the  right  hand  side  of  Equation  (2.55)  represents  residual  errors  after 
double  differencing.  For  baseline  separation  on  the  order  of  10  km  or  less,  double  difference 
residual  errors  are  on  the  order  of  centimeters  or  less  [37] . 
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2.8.3  Triple  Difference  Carrier-Phase  Differential  GPS.  The  previous  differenc¬ 
ing  process  could  be  extended  to  form  triple  differences.  These  triple  differences  take  the 
difference  between  at  time  tm,  and  at  some  later  time  Assuming 

continuous  carrier  tracking  is  maintained  on  the  interval  tm  <t<tn,  this  triple  difference 
would  eliminate  the  integer  ambiguity  altogether.  However,  the  triple  difference  measure¬ 
ment  doesn’t  carry  much  information  and  is  not  generally  used  to  obtain  a  differential 
carrier-phase  positioning  solution.  (It  is  sometimes  used  for  initialization  in  a  static  base¬ 
line  case.)  Because  triple  differencing  assumes  continuous  carrier  tracking,  it  can  be  useful 
for  detecting  breaks  in  carrier  tracking,  called  “cycle  slips”  [35]. 

2.9  Modified  Multiple  Model  Adaptive  Estimator  for  GPS  Carrier-Phase  Integer  Ambi¬ 
guity  Resolution 

The  technique  of  using  multiple  Kalman  filters  for  GPS  double  difference  carrier- 
phase  integer  ambiguity  resolution  was  introduced  early  in  the  GPS  carrier-phase  ambigu¬ 
ity  research  history  [21].  The  concept  involves  building  multiple  extended  Kalman  filters, 
each  based  on  a  given  candidate  integer  ambiguity  vector.  Then,  the  measurement  residual 
vectors  rj  and  residual  covariance  matrices  Aj  from  each  of  the  elemental  filters 
are  used  to  determine  which  elemental  filter  is  most  likely  to  contain  the  correct  ambiguity 
vector.  Two  researchers  [18,  19]  have  proposed  this  method  for  integer  ambiguity  reso¬ 
lution.  Their  research  is  summarized  next.  Recall  from  Section  2.2  that  using  extended 
Kalman  filters  in  an  MMAE  is  an  ad  hoc  extension  which  is  not  an  optimal  Bayesian 
estimator. 

The  first  proposal  for  using  multiple  Kalman  filters  for  GPS  carrier-phase  integer 
ambiguity  resolution  was  presented  in  1983  [19] ,  prior  to  GPS  becoming  fully  operational. 
It  proposed  the  use  of  an  MMAE®  to  increase  the  speed  of  the  carrier-phase  ambiguity 
resolution  process  in  static  survey  applications.  This  research  included  a  series  of  relatively 
simple  simulations.  Each  case  simulated  double  difference  measurements,  with  a  maximum 
of  4  GPS  satellites,  leading  to  a  maximum  of  3  carrier-phase  ambiguities  to  be  determined. 

®[19]  used  the  name  “Magill  Filter”  after  the  inventor  of  the  technique,  D.T.  Magill  [26] 
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In  the  most  demanding  case,  [19]  demonstrated  that  the  MMAE  could  be  used  to  determine 
GPS  carrier-phase  integer  ambiguities  in  a  fairly  simple,  yet  realistic,  three  dimensional 
GPS  simulation.  Since  search  region  limiting  techniques  had  not  yet  been  developed,  [19] 
relied  on  a  brute  force  search  of  all  possible  integer  ambiguities  within  5  integers  of  the 
correct  value  (a  total  of  11  possible  ambiguities  per  measurement),  which  was  arbitrarily 
set.  All  simulations  utilized  a  short  baseline  between  receivers,  and  ignored  many  common 
errors  including  atmospheric  effects  and  multipath.  The  results  were  encouraging,  showing 
a  simulated  ability  to  determine  the  integer  ambiguities  with  tens  of  minutes  of  measure¬ 
ment  data,  versus  approximately  an  hour  of  data  for  batch  processing  techniques  available 
at  that  time.  However,  due  to  the  brute  force  search  method,  [19]  correctly  concluded  that 
the  search  region  in  a  realistic  scenario  would  quickly  grow  too  large.  In  an  effort  to  reduce 
the  search  region,  [19]  derived  a  position  and  floating  point  ambiguity  estimating  Kalman 
filter  ^  to  provide  an  initial  ambiguity  estimate  to  the  MMAE.  The  use  of  a  floating  point 
filter  front  end  is  common  practice  today. 

This  initial  research  demonstrated  the  potential  to  use  the  MMAE  for  GPS  carrier- 
phase  integer  ambiguity  resolution.  However,  due  to  computing  power  available  at  the 
time  and  the  fact  that  GPS  was  in  its  infancy,  a  great  deal  of  work  remained.  It  was 
concerned  only  with  static  survey  applications,  and  did  not  address  the  issue  of  kinematics. 
Additionally,  it  used  only  fairly  simple  simulations  to  demonstrate  the  technique.  The 
technique  was  not  demonstrated  with  any  real  GPS  data. 


More  recently,  [17,  18]  proposed  a  modified  MMAE  for  kinematic  carrier-phase  ambi¬ 
guity  resolution.  In  [17,  18]  it  was  hypothesized  that  modeling  inaccuracies  would  cause  the 
measurement  residuals  to  be  large  for  all  the  MMAE  elemental  filters®.  More  specifically, 
recall  Equation  (2.41),  repeated  here  for  convenience: 

/ ,  N  /z(q)|a,z{q_i)  (zjja  =  afc,  Z  (tj_i)  =  Zj_i) pfc  (tj_i)  /r,  a-\\ 

Pk{ti)  =  -  (2.41) 

/z(ti)|a,Z(ii_i)  (^*1^  =  aj,  Z  (tj_i)  =  Zii-i)  Pj  (ti-i) 

1=1 


^The  floating  point  ambiguity  Kalman  filter  estimates  the  integer  ambiguity  values  without  constraining 
the  estimates  to  be  integers.  Therefore,  the  ambiguity  estimates  in  a  floating  point  Alter  are  not  typically 
integer  values. 

®This  problem  is  related  to  the  caution  mentioned  in  Section  2.3  about  adding  dynamics  pseudonoise 
during  Alter  tuning. 
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Equation  (2.41)  is  the  conditional  probability  density  function  for  the  parameter  a^,  given 
the  measurement  history  Zj_i.  Thus,  pk  (ti)  is  used  to  determine  which  carrier-phase 
integer  ambiguity  vector  (the  ambiguity  vector  is  the  parameter  in  this  application)  is 
correct.  The  conditional  densities  in  Equation  (2.41)  are  given  by  Equation  (2.47),  also 
repeated  here  for  convenience: 

/z(ti)|a,Z(ti_i)  (Cil^  Z  (tj_i)  Zj_i) 

(27r)  2  |Afc  (f.  )  |2  I  2  J 

where  (U)  Pfc  {t~ )  (ti)  +  Rfc  (ti)  is  the  covariance  matrix  of  the  residuals 

ffc  (t^)  =  Zj  —  H  (ti)  X  (t~),  ®  from  elemental  filter  k.  The  residual  covariance  A^  (t“)  is 
calculated  from  each  of  the  K  elemental  filters.  Equations  (2.41)  and  (2.47)  demonstrate 
that  the  MMAE  relies  upon  the  propagated  state  estimates  x  (t~ )  to  calculate  measure¬ 
ment  residuals  and  determine  which  ambiguity  vector  is  correct.  However,  the  authors 
in  [17,  18]  felt  the  state  dynamics  model  didn’t  have  enough  fidelity  to  enable  the  MMAE 
to  distinguish  between  residuals  from  the  elemental  filters.  If  the  dynamics  model  was 
poor,  compared  to  the  precise  measurements,  then  the  measurement  residuals  from  all 
elemental  filters  would  be  large,  causing  MMAE  failure  in  determining  the  correct  ambi¬ 
guity  set.  To  correct  this  potential  problem,  [17,  18]  proposed  using  “post-fit”  residuals  in 
the  hypothesis  test  portion  of  Figure  1.2.  These  “post- fit”  residuals  are  the  measurement 
residuals  calculated  after  measurement  incorporation  in  each  of  the  elemental  filters.  This 
modification  appeared  to  fix  the  perceived  modeling  noise  problem.  Table  2.5  summarizes 
the  changes  made  to  the  traditional  MMAE. 

The  modified  MMAE  of  Table  2.5  was  built  and  tested  using  kinematic  (i.e.,  moving) 
double  difference  carrier- phase  measurements.  These  tests  showed  the  modified  MMAE 
successfully  resolved  GPS  carrier-phase  integer  ambiguities.  The  primary  contribution  of 
this  dissertation  is  to  provide  a  mathematical  foundation  which  characterizes  and  general- 

®To  be  completely  correct,  [17,  18]  implemented  an  extended  Kalman  filter  with  measurement  residuals 
given  by  rk  (t~)  =  Zi  —  h  [x  (t~)  ,  ti) .  The  definition  above  is  used  to  be  consistent  with  the  derivation 
of  the  MMAE  based  on  a  linear  Kalman  filter,  as  presented  in  Section  2.3. 
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Table  2.5  MMAE  Modification  for  GPS  Carrier-Phase  Integer  Ambiguity  Resolution  [18 


Traditional  MMAE 

Modified  MMAE 

rfc  {t~ )  =  z  (ti)  -  h  [xfc  {t~ )  ,  ti] 

{tf )  =  z{ti)  -  h  [xfc  [tf )  ,  ti] 

Afc  {t~ )  =  Hfc  (ti)  Pfc  {t~ )  U.I  (ti)  -h  R  (ti) 

Afc  (t+)  =  Hfc  iU)  Pfc  {tf)  {ti)  +  R  {ti) 

izes  the  modified  MMAE.  While  the  modified  MMAE  was  designed  for  GPS  carrier-phase 
integer  ambiguity  resolution,  development  of  the  GRMMAE  is  not  restricted  to  any  partic¬ 
ular  application.  The  modified  MMAE  of  [17,  18]  is  one  specific  example  of  the  GRMMAE 
developed  in  Chapter  III. 

2.10  Chapter  Summary 

The  background  necessary  to  understand  the  GRMMAE  derivation  in  Chapter  III 
was  provided  in  this  chapter.  The  Kalman  filter  was  derived  as  the  basic  building  block 
for  the  Multiple  Model  Adaptive  Estimator.  The  extended  Kalman  filter  was  developed 
as  it  is  frequently  used  in  an  ad  hoc  extension  of  the  MMAE.  The  hypothesis  conditional 
probability  formula  for  the  MMAE  was  derived  and  several  MMAE  performance  enhance¬ 
ments  were  presented.  Since  this  dissertation  is  inspired  by  the  use  of  a  modified  MMAE 
for  GPS  carrier-phase  integer  ambiguity  resolution,  an  introduction  to  GPS  was  presented. 
Carrier-phase  integer  ambiguities  were  described  and  single  and  double  difference  carrier- 
phase  differential  GPS  were  discussed.  The  chapter  concluded  by  describing  the  modified 
MMAE  for  GPS  carrier-phase  integer  ambiguity  resolution.  The  contribution  of  this  dis¬ 
sertation  is  to  provide  a  sound  mathematical  formulation  which  generalizes  this  modified 
MMAE. 
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III.  Theory 


In  this  chapter,  the  Generalized  Residual  Multiple  Model  Adaptive  Estimator  is  de¬ 
rived.  A  “generalized  residual,”  which  is  a  linear  combination  of  traditional  and  post-fit 
residuals,  is  defined.  Properties  of  the  generalized  residuals  from  elemental  filters  based  on 
both  the  correct  and  incorrect  hypotheses  are  given  for  some  common  MMAE  applications. 
The  generalized  residuals  and  their  properties  are  then  used  to  re-derive  an  MMAE  based 
on  the  generalized  residual,  rather  than  the  traditional  residual.  The  newly  derived  MMAE 
is  called  the  Generalized  Residual  Multiple  Model  Adaptive  Estimator  (GRMMAE).  The 
GRMMAE  has  exactly  the  same  form  and  purpose  as  the  traditional  MMAE.  The  only 
functional  difference  between  the  two  is  the  formula  used  for  hypothesis  conditional  prob¬ 
ability  calculation.  This  derivation  assumes  the  GRMMAE  elemental  filters  are  linear 
Kalman  filters.  The  effect  of  using  an  extended  Kalman  filter  for  the  GPS  carrier  phase 
integer  ambiguity  resolution  application  is  discussed  in  Ghapter  IV. 

Each  section  in  this  chapter  provides  details  of  a  specific  research  contribution  re¬ 
sulting  from  this  dissertation.  The  entirety  of  this  chapter  is  the  primary  contribution  of 
providing  a  mathematical  foundation  and  generalization  for  the  modified  MMAE.  In  addi¬ 
tion  to  its  application  to  the  modified  MMAE,  the  derivation  provided  is  general  enough 
to  apply  to  a  traditional  MMAE  or  any  other  MMAE  using  a  generalized  residual. 

3.1  Generalized  Residual 

Begin  by  defining  the  generalized  residual  vector  as  a  linear  combination  of  traditional 
and  post-fit  residuals 


=  7r(t.  )  +  (1  -  y)r{tf) 

(3.1) 

where  r*  (ti)  is  the  generalized  residual,  r  is  the  traditional  Kalman  filter  residual, 
r  (tf^  is  the  post-fit  residual,  both  defined  in  Table  2.5,  and  7  is  a  real,  scalar,  user-defined 
GRMMAE  design  parameter.  Although  7  is  not  constrained  to  be  in  the  region  from  0  to  1, 
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inclusively,  7  G  [0, 1]  has  the  most  physical  meaning.  It  is  anticipated  that  most  engineering 
applications  will  use  7  G  [0, 1].  Finally,  throughout  this  dissertation,  a  superscript  >1=  denotes 
variables  based  on  the  generalized  residual.  The  *  notation  does  not  denote  the  complex 
conjugate  anywhere  in  this  dissertation. 

A  simple  derivation  yields  a  matrix  transformation  from  the  traditional  to  the  gen¬ 
eralized  residual 


=z{ti)  -  7H(ti)x(t“)  -  (1  -  7) 

=  z{ti)  -  7H(ti)x(t-) 

-  (1  -  7)H(tj)  x(t“)  +  K{ti)[z{ti)  -  H(tj)x(t-)] 
=  [!-(!-  ^)U{U)K{U)][z{U)  -  H(t,)x(t-)] 
r*  (ti)  =  T  (ti)  r  {t- ) 


(3.2) 


where  T  (tj)  =  [l  —  (1  —  7)  H  (tj)  K  (tj)]  is  the  transformation  that  interrelates  r*  (tj) 
and  r  (t~Y  The  transformation  T  (tj)  and  design  parameter  7  are  discussed  further  in 
Sections  3.4  and  3.7,  respectively. 


3.2  Generalized  Residual  Properties  for  the  Correet  Elemental  Filter 

Generalized  residual  vector  properties  are  derived  in  this  section,  assuming  that  the 
GRMMAE  elemental  filter  contains  the  correct  hypothesis.  From  Equation  (3.2), 

rl{ti)  =  Tk{ti)rk{t-) 

where  the  subscript  k  indicates  variables  associated  with  the  k^^  elemental  filter.  Since  the 
kf’^  filter  is  assumed  correct,  (t~)  is  Gaussian  [29].  The  generalized  residual  (ti)  is  a 
linear  function  of  (t~),  and  therefore,  also  Gaussian.  Thus,  (ti)  is  completely  defined 
by  its  mean  and  covariance.  Given  Equation  (3.2),  the  conditional  mean  of  the  generalized 
residual  from  the  correct  elemental  filter  is 
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=E[r*{ti)  I  a  =  afc,  Z{ti-i)  =  Zi_i] 

=  E  T{ti)r{t~)  I  a  =  afc,  Z{ti-i)  =  Zi_i 
=  Tk{ti)E  r{t~)  I  a  =  a^,  Z{ti-i)  =  Zi_i 
mr*(t,)  =  0  (3.3) 

where  the  final  line  above  comes  from  E  [r  (t“)  |  a  =  a^,  Z(fj_i)  =  Zj_i]  =  0  [29]. 

Notice  that,  similar  to  the  development  in  Section  2.3,  the  measurement  history  at  ti-i  is 
used  in  Equation  (3.3). 

Since  r|,  (U)  is  zero-mean,  the  conditional  covariance  is  given  by 

AliU)  =  E[r*{U)r*^  iU)  |  a  =  a^,  Z(t,_i)  =  Z,_i 

=  E  T  {ti)r  [t~)  {T  {ti)r  \  a  =  a^,  Z{ti-i)  =  Zi_i 

=  Tk{ti)E  r{t~)r'^{t~)  |  a  =  a^,  Z{ti-i)  =  Zi_i  (ti) 

=  Tfc  (ti)  Afc  (t^  )  (ti) 

AliU)  =[!-(!-  ^)Uk{ti)Kk{U)][Hk{U)Pk{t-)nl{ti)  +  K{ti)] 

[I  -  (1  -  (3.4) 

where  the  substitution 

E  r{t~)r'^{t~)  I  a  =  afc,  Z{ti-i)  =  Zi_i  =  [Uk  {ti)Pk  {t~)  ill  (U)  +  ^{U)] 

is  given  in  [30],  and  Tk{ti)  is  brought  outside  the  expectation  because  it  is  not  a  random 
matrix. 

Consider  the  following  specific  cases  for  Equation  (3.4) 
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A-liU) 


I  Hfc  (ti)  Pfc  (t .  )  (tj)  +  R(tj)  7  =  1 

\  [I  -  Hfc  (ti)  Kfc  (ti)]  [Hfc  (ti)  Pfc  (t-)  (ti)  +  R  (ti)]  [I  -  Hfc  (t,)  Kfc  {U)f  7  =  0 

(3.5) 

The  7  =  1  case  is  the  traditional  MMAE  and  the  result  of  Equation  (3.5)  exactly  matches 
the  traditional  filter-computed  residual  covariance  [30] 

Afc  (ti)  =  [Hfc  (ti)  Pfc  {t~ )  (ti)  +  R  (ti)] 

The  7  =  0  case  is  a  GRMMAE  using  post-fit  residuals.  This  result  is  easily  confirmed 
using  Equation  (3.2). 

3. 3  Generalized  Residual  Properties  for  Ineorreet  Elemental  Filters 

Generalized  residual  properties  for  the  correct  GRMMAE  elemental  filter  were  de¬ 
rived  in  the  previous  section.  A  similarly  general  derivation  cannot  be  accomplished  for 
the  incorrect  elemental  filters.  However,  for  specific  applications,  generalized  residual  prop¬ 
erties  from  the  incorrect  elemental  filters  can  be  derived.  This  is  accomplished  for  four 
different  types  of  MMAE  problems:  1)  a  measurement  bias  problem;  2)  an  unknown  equiv¬ 
alent  discrete-time  dynamics  noise  covariance  problem;  3)  an  unknown  dynamics  matrix 
problem;  4)  an  unknown  measurement  noise  covariance  problem.  The  procedures  shown 
in  these  examples  can  be  applied  in  a  similar  manner  to  other  types  of  problems  as  well. 

3.3.1  Measurement  Bias  Parameter.  When  the  parameter  to  be  estimated  is  a 
constant  measurement  bias  vector,  the  residual  properties  for  both  the  correct  and  incorrect 
elemental  filters  can  be  derived.  An  example  of  such  a  parameter  is  the  GPS  carrier- 
phase  integer  ambiguity  vector  (the  problem  originally  motivating  this  research).  Table 
3.1  summarizes  the  relationship  between  the  elemental  filters  for  the  key  matrices  and 
vectors  in  the  measurement  bias  scenario  under  the  previous  assumption  that  the  elemental 
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Table  3.1  Relationship  Among  Elemental  Filter  Matrices  and  Vectors  for  Unknown  Mea¬ 


surement  Bias  Application 


Same  in  All  Elemental  Filters 

Different  in  Each  Elemental  Filter 

z(tj) 

z  {ti) 

Bd  {ti) 

X  {t~)  (after  first  measurement) 

u(ti) 

Qd  (ti) 

r*  {ti) 

H(t,) 

R(tj) 

K{U) 

A*  {ti) 

T(t*) 

$  {ti,ti-i) 

filters  are  based  on  linear  system  models.  Chapter  IV  discusses  the  impact  of  using  a 
nonlinear  measurement  model  for  GPS  carrier-phase  integer  ambiguity  resolution  on  the 
(linear  model)  derivation  presented  here. 

Assume  the  correct  parameter  is  contained  in  the  elemental  filter.  The  sensor 
measurement  can  be  modeled  by 


z  (ti)  =  Zd  (ti)  +  rifc  (3.6) 

where  z  (ti)  is  the  measurement  output  from  the  sensor,  Zd  (ti)  is  the  desired  (unbiased) 
measurement,  and  is  the  unknown,  constant  measurement  bias  vector  (i.e.,  the  param¬ 
eter  to  be  determined).  Note  that  the  desired  measurement  Zd{ti)  is  a  noise-corrupted 
measurement.  It  is  termed  “desired”  here  because  it  is  unbiased. 

For  ease  of  implementation,  each  elemental  filter  can  be  designed  to  operate  on  a 
corrected  measurement  given  by  [18,  20] 


2jr  (^j)  —  2  (^i)  rijr 

Zj{ti)  =  Zd{ti)  -F  rifc  -  THj  (3.7) 
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where  rij  is  the  constant  measurement  bias  hypothesized  in  the  elemental  filter  and 
is  assumed  to  be  the  correct  bias. 

Using  the  corrected  measurement  of  Equation  (3.7)  leads  to  all  elemental  filters  hav¬ 
ing  the  same  structure  but  operating  on  different  measurements.  It  is  particularly  impor¬ 
tant  to  note  that 


Hfc  {ti)  =  Hj  {ti 

)  =  H(U 

and 

Kfc  {ti)  =Kj  {ti 

)  =  K(U 

(3.8) 


Equation  (3.8)  is  valid  because  of  the  way  the  elemental  filters  were  constructed  for 
this  application.  It  implies  that  the  transformation  matrix  T  (U)  is  the  same  in  every 
elemental  filter. 

Assuming  the  measurement  history  up  to  and  including  U_i  is  available,  the  residual 
vector  for  any  elemental  filter,  given  Equation  (3.7),  is 


Y*{ti)  =Zj{ti)  -  7H(ti)%(t“)  -  (1  -  7)H(U)%(t+) 

=  Zd{ti)  -f  rifc  -  rij  -  7H(ti)xj(t“) 

-  (1  -  7)H(U)|%(t“)  -F  K(U)[zrf(U)  -f  rifc  -  rij  -  H(U)%(t“)]| 
r*  (U)  =  Zd  (ti)  -f  Arifcj  -  7H  (ti)  Xj  {t~ ) 

-  (1  -  7)  H  (U)  {Xj  (t“)  +  K{ti)[zd{ti)  +  Aukj  -  H(U)Sj(t-)]| 

(3.9) 

where  An^j  =  —  rij  is  the  difference  between  the  hypothesized  measurement  bias  in 

elemental  filter  j  and  the  correct  measurement  bias  in  elemental  filter  k.  Equation  (3.9) 
can  be  simplified  to 
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rj  (ti) 


r*j  {ti) 


=  [I  -  (1  -  7)H(t,)K(t,)][zrf(tO  -  +  Arifc,] 

=  T  (ti)  [zd  {ti)  -  H  {ti)  Xfc  {t~ )  +  H  {ti)  Xfc  {t~ )  -  H  {ti)  Xj  {t~ )  +  Arifcj] 

=  T{ti)[rk{t-)  +  H(ti)(Xfc(t“)  -  ^j{t~))  +  Arifcj] 

_  f  T  {ti)  Yk  {t~ )  j  =  k 

\  T  {ti)  [xk  {t- )  +  H  {ti)  (Xfe  {t- )  -  %  {t- ))  +  Arifcj]  Vj 

(3.10) 


In  the  third  line  of  Equation  (3.10),  is  equated  to  Zd{ti)  —  H(tj)xfc  (t~). 

That  relationship  only  holds  for  the  (i-e.,  correct)  elemental  filter.  In  general,  Yj  (t“)  = 
Zj  {ti)  —  H  {ti)  Xj  {t~).  In  the  special  case  when  j  =  k,  zj  {ti)  =  Zd  {ti),  as  seen  in  Equation 
3.7,  which  leads  to  {t~ )  =  Zd  {U)  -  H  {ti)  x*.  (t~ ) . 

Given  E  {t~)  \  a  =  a^,  Z(tj_i)  =  Zi_i  =  0  [29]  and  Equation  (3.10),  the 
expected  value  of  r*  {ti)  is 


nir*{ti)  =E\Y*{ti)  I  a  =  afc,  Z{ti-i)  =  Zi_i 

f  0  j  =  k 

mr*  {ti)  =  I  (3.11) 

[  T  {ti)  [H  {ti)  (Xfc  (f.  )  -  Stj  (tj  ))  +  Arifcj]  Vj 

The  covariance  of  r*  (ti)  is 


A*  (ti)  =E  Y*{ti)Yf  {ti)  I  a  =  afc,  Z(ti_i)  =  Zi_i  -  rrir*  (t*)  (ti) 

=  T  (ti)  e|  [xk  (t“ )  +  H  (ti)  (Xfc  (t“ )  -  Xj  (t"))  +  Arifcj] 

[rA:  (t“) +H(ti)  (xfc  (t“) -Xj  (t“))  +  Arifcj]^  I  a  =  a^,  Z(ti-i)  =  Zi_i|T'^(ti) 

-  mr*  (ti)  (ti) 

^  J 
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A*  (ti)  =  T  (ti)  Yk  (tj  )rl  {t-)  I  a  =  a^,  Z  (ti-i)  =  Zi_i 


+ 

E 

I'fc 

{t'i 

-)  1  a  = 

Z  {ti— 

i)  =  2 

'h-1 

(xfc  {ti  )  {ti  ))'^H(t, 

■f 

+ 

E 

I'fc 

-)  1  a  = 

2  {ti— 

i)  =  2 

'h-i 

Ani; 

+ 

H 

{ti 

)  (Xfc  {ti  )  - 

%  {ti))E 

•■fc  {ti ) 

a 

—  ^ki  Z  {ti—l)  —  1 

+ 

H 

{ti 

)  (Xfc  )  - 

%  {ti  ))  (xfc  {ti  )  - 

-% 

{ti)y 

{ti)^ 

+ 

H 

{ti 

)  (Xfc  [tf )  - 

%  {h  ))A: 

T 

^kj 

+ 

AukjE 

{ti) 

3  = 

Z(t*-i) 

= 

Zi-i 

+  ^ 

An 

kj  ( 

Xfc 

{ti )  -  ^ 

{U)^  + 

Aukj An^j'^T  {ti)  -  inr*{t, 

i)  mf*  {ti) 

3 

(3.12) 

The  mean  and  covariance  of  the  traditional  residual  from  the  correct  elemental  filter 
Tfc  (t“)  are  used  to  simplify  Equation  (3.12): 


A*  (ti)  =  T  (ti)  |Afe  {t~ )  +  [H  (ti)  (Sfc  {t~ )  -  Xj  {t~))  +  Aukj] 

[H(tj)  (Sfc  (t“)  (t“))  +  Anfcj]^}T^(fi)  - 

=  Afc  {ti)  +  T  {ti)  [H  {ti)  (Xfc  (t“ )  -  %  (f“))  +  Aukj] 

[H  {ti)  (Xfc  (f“ )  -  Sj  (t“ ))  +  Aukj]  (tj) 

-  T  (ti)  [H  (ti)  (Sfc  (t“ )  -  %  (t“ ))  +  Arifcj] 

[H  (ti)  (Xfc  (t“ )  -  %  (t“ ))  +  Arifcj]  (ti) 

A*(t,)  =A^(t,)  =  T(t,)Afc(t-)T^(t,)  (3.13) 


Thus,  the  generalized  residuals  from  all  elemental  filters  will  have  the  same  filter- 
computed  covariance  and  differing  means  when  the  parameter  of  interest  is  a  measurement 
bias. 


3.3.2  Residual  Properties  for  Varying  Qd{ti)  Parameter.  Similar  to  the  previ¬ 
ous  subsection,  residual  properties  can  also  be  established  in  the  more  traditional  MMAE 
application  of  determining  the  correct  equivalent  discrete-time  dynamics  driving  noise  co- 
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Table  3.2  Relationship  Among  Elemental  Filter  Matrices  and  Vectors  for  Unknown 


Application 


Same  in  All  Elemental  Filters 

Different  in  Each  Elemental  Filter 

^  {ti,ti-i) 

X  {tJ ) 

Bd  {ti) 

^W) 

u(ti) 

H(t,) 

R(tj) 

K(t,) 

z(tj) 

r*  {ti) 

A*  {ti) 

TiU) 

Qd  {ti) 

variance  (tj).  In  this  problem,  each  elemental  filter  is  built  with  a  different  hypothesized 
Qrf  (tj)  matrix.  Table  3.2  summarizes  the  relationship  among  the  key  elemental  filter  ma¬ 
trices  and  vectors  for  this  scenario,  assuming  that  the  elemental  filter  models  are  linear. 

Assume  elemental  filter  k  is  based  on  the  correct  hypothesis.  To  calculate  the  mean 
of  r*  (tj)  first  define  r*  (tj)  in  terms  of  {t~) 


r*  {k 

)  =  [I  - 

(1  -  7)H,(tOK, 

{ti)h  ih ) 

= 

0  U  {ii ) 

r*  {ti 

)  = 

0  {ti)  +  H  {ti) 

(Sfc(f-)  - 

(3.14) 

Given  Equation  (3.14),  the  mean  and  covariance  of  r*  (fj)  are  calculated  from  the 
mean  and  covariance  of  (t~ ) . 


m,*  {ti) 

m,*  {ti) 

A*  (ti) 

A*  iU) 


E 

E 


Y*{ti)  I  a  =  afc,  Z{ti-i)  =  Zi_] 


Tj(t*)  Yk{t^)  +  H(ti)(xfc(t.)  -  Xj(tj))  I  a  =  afc,  Z{ti-i) 


—  Tj  {ti)  H  {ti)  Xfc  (t^  )  Xj  {t^  ) 


=  E 


r*  {ti)  r*/  {ti)  I  a  =  afc,  Z  {U-i)  =  Zi_i  -  rrir*  {U)  {U) 

■I  I  I  3 


—  Tj  {ti)  Afc  {t^  )  T j  {ti) 


(3.15) 
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As  in  the  measurement  bias  application,  the  generalized  residual  vectors  will  have 
differing  means,  dependent  upon  the  elemental  filter.  However,  unlike  the  measurement 
bias  application,  the  filter-computed  residual  covariance  will  also  vary  among  the  elemental 
filters  because  the  transformation  matrix  (tj)  varies  among  the  elemental  filters. 

3.3.3  Residual  Properties  for  Varying  F  (t)  Parameter.  The  linear  difference 
equation  for  system  dynamics  in  Equation  (2.1)  is  calculated  from  the  differential  equation 
[29] 


x(t)  =F(t)x(t)  -|-  B(t)u(t)  -|-  G  (t)  w  (t)  (3.16) 


where  F  (t)  is  the  system  dynamics  matrix.  The  state  transition  matrix,  $(tj,tj_i),  in 
Equation  (2.1)  is  the  matrix  which  satisfies  the  differential  equation  and  initial  condition: 


d[^  {t,to)] 
dt 


=  F  (t)  $  {t,to) 


^  {to,  to)  =  I 


(3.17) 


For  a  time  invariant  system  (i.e.,  constant  F  (t)),  the  state  transition  matrix  can  be  calcu¬ 
lated  using  the  matrix  exponential 


^{t,to)  =  ^{t-to)  =  (3.18) 

Thus,  the  system  dynamics  matrix  F  {t)  is  used  to  define  the  linear  difference  equation  of 
Equation  (2.1).  With  this  in  mind,  an  example  of  a  GRMMAE  with  an  unknown  F  {t)  is 
presented  next. 

In  this  scenario,  it  is  assumed  that  the  matrix  describing  the  system  dynamics,  F  {t) 
contains  an  unknown  parameter.  Such  a  case  arises  when  the  F  (t)  matrix  contains  a 
time  constant,  and  this  approach  useful  for  target  identification.  Different-valued  time 
constants  are  hypothesized  for  various  targets,  with  an  observed  target  being  classified  by 
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Table  3.3  Relationship  Among  Elemental  Filter  Matrices  and  Vectors  for  Unknown  F  (tj) 


Application 


Same  in  All  Elemental  Filters 

Different  in  Each  Elemental  Filter 

u(ti) 

Bd  {U) 

H(ti) 

R(ti) 

^  {ti,ti-i) 

^{ti) 

Qd  {ti) 

K(t,) 

^W) 

TTfl 

T(t*) 

r*  {ti) 

determining  the  most  probable  time  constant.  Since  the  state  transition  matrix  $  (U,  ti-i), 
is  calculated  directly  from  F  (t),  this  scenario  is  equivalent  to  an  unknown  $  (U,  t)  scenario. 
Additionally,  since  (tj)  is  calculated  from  $  (tj,  r)  using 


U 

Qrf(ti_i)  =  j  $(t*,r)G(r)Q(r)G'^(r)$^(ti,r)(iT  (3.19) 

ti-i 

this  scenario  is  related  to  the  unknown  (tj)  scenario  of  Section  3.3.2.  Table  3.3  lists  the 
relationship  of  the  key  elemental  filter  matrices  and  vectors  for  this  scenario.  Under  the 
conditions  in  Table  3.3,  specifically  each  filter  operating  on  the  same  measurement  vector 
and  the  transformation  Tj  (tj)  varying  between  the  elemental  filters,  the  derivation  of  the 
mean  and  covariance  for  r*  {ti)  is  exactly  the  same  as  in  Section  3.3.2.  Thus,  Equation 
3.15  applies  to  this  scenario. 

3.3.4  Residual  Properties  for  Varying  R(ti)  Parameter.  A  third  parameter  es- 
timatable  using  an  MMAE  is  the  measurement  noise  covariance  matrix  R(U).  In  this 
problem,  multiple  measurement  noise  covariance  matrices  are  hypothesized  in  the  elemen¬ 
tal  filters.  The  MMAE  determines  which  elemental  filter  contains  the  correct  R  (U)  matrix. 
Such  an  application  arises  in  failure  detection  if  the  noise  in  a  sensor  changes  as  the  result 
of  a  failure  (or  partial  failure).  This  application  also  arises  when  hypothesizing  differ- 
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Table  3.4  Relationship  Among  Elemental  Filter  Matrices  and  Vectors  for  Unknown 


R  (ti)  Application 


Same  in  All  Elemental  Filters 

Different  in  Each  Elemental  Filter 

$  (tj,tj_i) 

X  {t- ) 

Bd  (tj) 

u(ti) 

H(tj) 

TTfl 

Qd  (ti) 

K(t,) 

ziU) 

r*  (ti) 

A*  (tj) 

TiU) 

R(tj) 

ent  levels  of  interference  or  jamming  in  GPS.  Table  3.4  shows  the  relationship  of  the  key 
Kalman  filter  parameters  among  the  elemental  filters.  Under  the  conditions  in  Table  3.4, 
specifically  each  filter  operating  on  the  same  measurement  vector  and  the  transformation 
matrix  Tj  (ti)  varying  between  the  elemental  filters,  the  derivation  of  the  mean  and  co- 
variance  for  r*  (ti)  is  exactly  the  same  as  in  Section  3.3.2  and  Equation  (3.15)  applies  to 
this  scenario  as  well. 

3-4  Effect  of  T  (ti)  Transformation 

The  transformation  T  (t*)  =  [l  —  (1  —  7)  H  (tj)  K  (tj)]  occurs  frequently  in 

derivations  using  r*  (tj).  To  understand  the  effect  of  T  (tj)  on  the  traditional  residual 
r  (t~),  define  the  difference  Ar  (tj)  as: 


Ar  (tj)  =  r  (t .  )  -  r*  (tj)  (3.20) 

Using  Equation  (3.2),  it  can  be  shown  that 


Ar(tj)  =  (1  -  7)  H  (tj)  K  (tj)  r  (t- ) 


(3.21) 
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Examination  of  the  Kalman  filter  state  update  equation  from  Table  2.1,  repeated 


here  for  convenience, 


+  K(ti)[zi  -  H(ti)S(tJ] 

reveals  that  the  Kalman  gain  matrix  K  (tj)  transforms  measurement  domain  vectors  into 
the  state  space.  It  also  shows  that  the  measurement  matrix  H  (tj)  transforms  state  vec¬ 
tors  into  the  measurement  domain.  Thus,  multiplication  by  H  (tj)  K  (tj)  transforms  the 
measurement  domain  residual  vector  into  the  state  space,  then  back  into  the  measurement 
space.  A  two  dimensional  representation  is  shown  in  Figure  3.1.  The  scalar  (1  —  7)  is 
thought  of  as  a  gain  which  determines  how  much  the  traditional  residual  vector  is  scaled. 
After  scaling  by  (1  —  7),  the  residual  vector  is  rotated  (and  possibly  scaled  further)  by  the 
H  {ti)  K  (ti)  multiplication.  In  the  extreme  case  of  7  =  1  (i.e.,  r*  (U)  =  r  the  scalar 

gain  is  0,  and  Ar  (ti)  is  0.  In  general,  T  (ti)  is  a  transformation,  not  a  projection,  since 
T  (ti)  is  not  idempotent.  This  fact  was  confirmed  using  a  simulation  to  calculate  T  (ti)  in 
a  measurement  bias  application.  Although  a  formal  proof  that  T  (ti)  is  never  a  projection 
was  not  pursued,  it  is  expected  that  T  (ti)  will  not  be  a  projection  except  possibly  in  very 
unusual  cases. 

3.5  Hypothesis  Conditional  Probability  Derivation  Using  Generalized  Residuals 

This  section  is  similar  to  Section  2.3,  deriving  the  probability  that  each  GRMMAE 
elemental  filter  is  correct,  based  on  a  linear  combination  of  the  traditional  Kalman  filter 
residual  and  the  post-fit  residual.  To  calculate  the  probability  that  a  given  elemental  filter 
contains  the  correct  parameter  a,  consider  the  conditional  probability  density  function 

/a|r*(7),  Z(7_i)  (ti)  =  p*i,Z  (U-i)  =  Zi_i)  (3.22) 

This  is  equivalent  in  information  content  to  /z(ti)|a,z(ti_i)  (CJ^  =  a,  Z  (ti-i)  =  Zj_i)  from 
Equation  (2.37).  Equation  (3.22)  is  used  to  derive  a  hypothesis  conditional  probability  for¬ 
mula  similar  to  Equation  (2.41),  but  acting  upon  the  generalized  residual,  rather  than  the 
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(1-7)K(«,) 


Figure  3.1  2-D  Representation  of  Difference  Between  Traditional  and  Generalized  Resid¬ 

ual  Vectors 

traditional  residual.  Similar  to  Section  2.3,  assume  a  can  take  any  value  in  the  continuous 
range  A  C  M^.  Applying  Baye’s  rule  to  Equation  (3.22)  yields 


/a|r*(q),  Z(q_i)  Hr*  {U)  =  p*,Z  {U-i)  =  Zi_i)  = 

fr*{ti)\a,  Z(ii_i)  (P  1^  ^  (R— l)  Zj_i)  /a|Z(ti_i)  (^i— l)  Zj_i) 

//r*(q)|a,  Z{ti_i)  (pla  =  cx,Z{ti-i)  =  Zi_i)/a|z(q_i)  HZ{ti-i)  =  Zi-i)da 
A 

(3.23) 

Next,  for  feasible  implementation,  assume  the  parameter  vector  is  restricted  to  be 
a  member  of  the  discrete  set  {ai,  a2,  ...  ,ai^}.  Additionally,  assume  that  the  correct 
parameter  value  is  an  element  of  that  discrete  set.  Then,  define  the  hypothesis  conditional 
probability  as 


PkiU)  =  prob{a  =  ak\r*  {ti)  =  r*f.  {ti) ,  Z  (tj_i)  =  Zi_i}  (3.24) 
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Under  these  conditions,  the  conditional  hypothesis  probability  can  be  calculated  by 


evaluating  Equation  (3.23)  using  the  realization  of  the  measurement  history  Zj_] 


Pkiti)  = 


/r*(q)|a,  Z(ti_i)  (U)  1^  3-A:)  Z  Zj_i) 

K 

//r*(q)|a,  Z(q_i)  (r*  (U)  |a  =  «,Z(U_i)  =  Zi_i)  X]  Pa|Zi_i<^  (a  “  a^-)  da 

A  J  =  1 

/r*(q)|a,  Z(ti_i)  (^*)  1^  ~  ^  ~  Zj-i) 


A 


E  /  /r*(id|a,  Z(q_i)  (r*  (U)  |a  =  Q,  Z  (U_i)  =  Zi_i)  Pa|Zi_i<^  (« 

3  =  lA 

/r*(ii)|a,  Z(q_i)  {^k  {ti)  |3  =  afc,  Z  (U-l)  =  Zi_i)  Paj^|Zi_i 
E  /r*(q)|a,  Z(ti_i)  fl'j  (u)  |a  =  Uj,  Z  (U-l)  =  Zj_i'j  Paj|Zi_i 

j  =  1  ^  / 


aj)  da 


(3.25) 


where  Pafc|Zi_i  is  the  discrete  probability  mass  associated  with  the  parameter  realization 
given  the  measurement  history  realization  Zj_i.  To  see  that  Pafc|Zi_i  =  Pfc  (U-i),  consider 


~  ^*afc|zi_i-[I  -  (1  -  7)Hfc(ti)Kfc(q)]rfc(q"_j),  Zi_2 

P^k\^'l,{^i  —  l)^  Zi_2 

Pafc|Zi_i  =Pfc(t*-l)  (3-26) 

Substituting  Equation  (3.26)  into  Equation  (3.25)  yields  the  final  recursive  relation¬ 
ship 


Pkiti)  = 


fr*(u)\a,  Z(q_i)(rfc  iu)  |a  =  afc,Z(ti_i)  =  Z^_i)p*  (U_i) 

E  /r*bi)|a,  Z(ti_i)  (r^  (ti)  |a  =  aj,  Z  (U_i)  =  Zi_i)  p*  (U_i) 
i  =  1  ^  ^ 


(3.27) 


In  Equation  (3.27),  the  density  A* (t.)|a,  z(ti_i)  {p*  iU)  |a  =  a,  Z(U_i)  =  Zi_i)  is  Gaus¬ 
sian.  To  demonstrate  this,  consider  fz{ti)\a,  z(ti_i)  (Ci|a  =  a,  Z(U_i)  =  .Zj_i).  Erom  Sec¬ 
tion  2.3  and  [30],  this  second  density  is  Gaussian.  As  shown  in  Equation  (3.2),  the  gen- 
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eralized  residual  r*  (ti)  is  a  linear  function  of  z(fj)  and  x(f“).  When  Equation  (3.1)  is 
conditioned  on  Z(tj_i)  =  Zj_i,  as  in  Equation  (3.27),  x(t“)  becomes  the  realization 
x(t“).  Under  this  condition,  r*  (ti)  is  a  linear  function  of  the  random  variable  z(tj), 
and  the  realization  x  {t~).  Therefore,  since  /z(q)|a,  z(q_i)  (Cl^  =  ct,  Z(tj_i)  =  Zi^i)  is 
Gaussian  [30],  /r*(q)|a,  z{ti_i)  (Pfcl^  =  a,  Z  (tj-i)  =  .Zj_i)  is  also  Gaussian.  The  mean  and 
covariance  of  r*  (ti),  derived  in  Section  3.1,  completely  specify  this  conditional  density. 

3.6  Relationship  Between  pk  (U)  andp\{ti) 

The  traditional  MMAE  (7  =  1)  calculates  elemental  filter  conditional  probabilities 
from  Equation  (2.41),  repeated  here  for  convenience: 

Pk{ti)  =  -  (2.41) 

S  fz(ti)\a,  Z(ti_i)  (^*1^  =  Z  (tj-i)  =  Zj_i)  Pj  (tj_i) 

j  =  1 

Equation  (3.27)  and  Equation  (2.41)  must  evaluate  to  the  same  value  when  7  =  1 
(i.e.,  the  generalized  residual  becomes  the  traditional  residual).  Therefore,  there  is  a 
relationship  between  these  two  equations. 

Section  3.2  shows  that  (ti)  is  Gaussian  and  zero-mean,  with  covariance  given  by 
Equation  (3.4)  (assuming  the  elemental  filter  contains  the  correct  parameter).  Thus, 
the  conditional  density  is 


frl{ti)\a,  Z(q_i)  i^k  i^i)  1^  ~  ^  (^*-l)  ~  Zj_i) 


I  -1 


=  PUU)eM-2^l  it^)[^k{U)A^^{t-)Tl{U)\  rl{U)}  (3.28) 


In  Equation  (3.28),  (ti),  Ak  (tj  )  and  T^  (ti)  are  given  by 
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1 


Afc  (t- ) 
Tfc  {ti) 


{27r)f\Tk{U)[llk{ti)Pk{t-)lll{ti)  +  Kk{U)]Tl{ti)\'2 
[Hfc  (ti)  Pfc  {t~ )  Ul  (ti)  +  Rfc  {ti)] 

[I  -  (1  -  7)  Hfc  (ti)  Kfc  (ti)]  (3.29) 


Consider  the  term  inside  the  exponential  first 


-Jrf  (i*)Afc  ^(ti)4(ti) 


1 
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=  --4  (t,)  Tk{U)A-^t-)Tl{U)  rliU) 


~  2^^  (^i)]'^fc  i^i)  (^*)  {ti  ) 

=  -  ^rl  {t-)Tl  {U)[Tl  {U)]-\Uk{U)Pk{t-)Ul  iU)  +  Rfc(tO]“'T-i(t,)Tfc(t,)rfc(t-) 
1 


=  -  x^fc  (^i  )  [Hfc  (ti)  Pfc  (t.  )  (ti)  +  Rfc  (ti)]  rfc  (t.  ) 


=  -  ih  )  Afc  ^  {ti  )  rfc  (t,  ) 


(3.30) 


Comparing  Equation  (3.30)  to  Equation  (2.47)  shows  that  the  term  inside  the  expo¬ 
nential  for  /r*(i.)|a,  z(7_i)  (rfc  {ti)  |a  =  afc,  Z  (ti_i)  =  Zi_i)  is  exactly  the  same  as  the  term 
inside  the  exponential  for  /z(i7|a,  z(7_i)  (zi|a  =  afc,Z(ti_i)  =  Zi_i).  This  implies  that  the 
GRMMAE  design  parameter  7  has  no  effect  on  the  likelihood  quotients.  Thus,  regardless 
of  the  value  of  7,  the  likelihood  quotient  inside  the  Gaussian  densities  of  Equation  (3.27) 
is  —  l/2r^  {ti)  {ti)  rfc  (tj),  the  same  as  it  is  in  the  traditional  MMAE. 

Now,  consider  the  {ti)  term  in  front  of  the  exponential  in  Equation  (3.28).  The 
determinant  |Tfc  (tj)  [Hfc  (tj)  Pfc  (t~)  H^  {ti)  +  R(ti)]T^  {ti)  \  can  be  broken  into  the  3 
separate  determinants  in  the  first  line  of  Equation  (3.31). 


/3fc  {ti) 


1 

(27r)^{|Tfc(t,)|  |[Hfc(t,)Pfc(t-)Hnt*)  +  Rfc(t*)]|  \Tl{U)\y^ 
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h’k  - r 

(27r)^{|Tfc(ti)|  I  [Hfc  (to  Pfc  (t“)  +  Rfc(t*)]|  |Tfc(tOI}" 

1 

(27r)?abs(|Tfc(tOI){  |  [H^  (tO  Pfc  (t“)  HHi*)  +  Rfc  (ti)]  I  }" 

(51{U)  =abs(|T,(tOI)“' - — - 

(27r)^{  \[Uk{U)Pk{t-)lll{U)  +  Rfc(tO]| 


(3.31) 


In  Equation  (3.31)  the  determinant  jT^  (tj)  |  is  squared  inside  the  square  root,  then  brought 
outside  the  square  root.  When  it  is  brought  outside,  the  positive  square  root  is  used. 
Therefore,  the  notation  abs(|Tfc  {U)  |)  is  used  to  denote  a  positive  square  root. 

Equation  (3.31)  equals  the  beta  term  from  /z(q)|a,  Zi_i  (zi|a  =  afc,Z(tj_i)  =  Zi_i), 
in  Equation  (2.47),  multiplied  by  the  scalar  term  abs(|Tfc  (t*)  |~^).  Thus, 


frl(ti)\a,  Z(ti_i)  {'^k  i^i)  1^  ~  Z  (tj-l)  —  Zj_i) 

=  abs(|Tfc(tO  l)^Vz(q)|a,  Z{q_i)  (z(ii)  |a  =  Z  (ti_i)  =  Zi_i) 

(3.32) 


Under  the  condition  7  =  1  (i.e.,  (U)  =  Vk  )),  T^  (U)  =  I  and 


frl{ti)\a,Z{ti-i)  (^k  i^i)  1^  ~  afc,  Z  (U-l)  —  Zi_i) 

=  /z(q)|a,z{ti_i)  (z(U)  |a  =  afc,Z(U_i)  =  Zi_i)  (3.33) 

Thus,  when  7  =  1,  the  GRMMAE  is  equivalent  to  a  traditional  MMAE. 

Equations  (3.30)  through  (3.32)  demonstrate  that  the  GRMMAE  design  parameter 
7  affects  only  the  /3*  (U)  normalization  term  in  front  of  the  exponential  in  the  Gaussian 
densities  of  Equation  (3.27).  Thus,  7  in  the  generalized  residual  has  the  effect  of  artificially 
increasing  or  decreasing  the  (5*  {U)  term  magnitude  in  these  Gaussian  densities.  Ghapter 
IV  demonstrates  that  this  results  in  a  susceptibility  to  the  /3-dominance  effect. 
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3.1  Optimization  of  p\  {ti) 

The  objective  of  the  7  parameter  in  the  generalized  residual  is  to  improve  the  rate  of 
convergence  to  the  correct  elemental  filter  of  the  GRMMAE  compared  to  the  traditional 
MMAE.  Thus,  under  the  assumption  that  the  elemental  filter  is  correct,  the  optimal 
value  of  7  should  maximize  (tj)  given  in  Equation  (3.27).  All  GRMMAE  applications 
can  be  divided  into  2  groups  to  accomplish  this  objective.  The  first  group  consists  of  ap¬ 
plications  in  which  the  transformation  matrix  T  (tj)  is  identical  in  all  elemental  filters.  Eor 
this  scenario,  any  7  may  be  used  with  identical  results  achieved.  This  will  be  demonstrated 
in  Section  3.7.1.  The  second  group  consists  of  applications  in  which  the  transformation 
matrix  T  (tj)  varies  among  the  elemental  filters.  In  these  applications,  a  specific  7  value 
will  maximize  (t*).  The  specific  optimal  value  of  7  will  vary  dependent  upon  the  ap¬ 
plication.  However,  these  applications  will  be  susceptible  to  /3-dominance  if  7  /  1,  as 
demonstrated  in  Sections  3.7.2  and  4.3. 

In  both  cases,  the  same  general  procedure  is  used  to  determine  the  optimal  7  value: 

1.  Determine  the  generalized  residual  t*-  (tj)  in  terms  of  the  traditional  residual  for  the 
correct  elemental  filter  {t~). 

2.  Determine  the  covariance  of  the  generalized  residual  A*  (tj)  in  terms  of  the  covariance 
of  the  traditional  residual  for  the  correct  elemental  filter  {t~). 

3.  Perform  algebraic  substitution  and  simplification  to  express  the  generalized  hypoth¬ 
esis  conditional  probability  of  Equation  (3.27)  in  terms  of  (t“)  and  A^  (t~). 

4.  Determine  the  7  value  which  maximizes  p^  (tj)  (either  in  steady  state  or  at  each 
measurement  update). 

Steps  1  and  2  have  already  been  accomplished  in  Sections  3.2  and  3.3.  Because  of  the 
/3-dominance  effect.  Step  4  will  not  be  pursued.  Rather,  the  designer  is  strongly  cautioned 
against  using  7  7^  1  unless  T  (tj)  is  identical  in  all  elemental  filters.  The  reason  for  this 
caution  will  be  demonstrated  in  Sections  3.7.2  and  4.3. 

3.7.1  Same  T  (U)  In  All  Elemental  Filters.  In  this  section.  Equation  (3.27)  is 
applied  to  applications  with  T  (t*)  identical  in  all  elemental  filters.  Significantly  general 
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notation  and  equations  applicable  to  all  cases  with  identical  T  (tj)  in  all  elemental  filters 
cannot  be  developed.  Therefore,  the  specific  example  of  estimating  a  constant  measurement 
bias,  as  detailed  in  Section  3.3.1,  is  shown.  This  will  demonstrate  the  4-step  procedure 
outlined  in  the  previous  section.  Other  applications  with  the  same  T  (tj)  in  all  elemental 
filters  can  be  similarly  derived.  Begin  with  Equation  (3.27),  substituting  the  residuals  from 
the  elemental  filters  into  the  densities 

Pk  iU)  =  - ^ ^ - 


(3.34) 


where  the  numerator  comes  from  Equation  (3.32) 

Next,  use  Equation  (3.13)  to  substitute  for  A*  (t*) 


Pkiti)  = 


(27r)' 


•'%bs(|T|‘..)|)|A.(,. -)!■/»  «p{-i4(‘r)A,-‘(i-)r^(t-)} 


K 

y  _ 1 _ 

(2^)"*/^abs(|T(q)|)|Afc(q-)|i/2 


exp 


{-irf  iU)  [T  iU)  Afc  (t-)  T^  iU)]  r*  {U)}  p*  (t._i) 


(3.35) 


Notice  in  Equation  (3.35)  that  all  densities  in  the  both  the  numerator  and  denom¬ 
inator  are  multiplied  by  1/  (27r)™''^^  abs(|T  (t*)  |)  |  A^  (f“)  .  Therefore,  this  term  is 

eliminated 


Pk  (ti) 


Eexp{g*}p*(ti_i) 
1=1  ^  ^ 


(3.36) 


where 
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Qj  =  -  \  [rfc  {t~ )  +  H  (ti)  {xk  {t- )  -  Xj  (t“))  +  Arifcjj'^T^  (ti)  A*  '  (U)  T  (ti) 

[vk  {t- )  +  H  (ti)  (Sfc  {f- )  -  %  (t-))  +  Aukj] 

=  {t~ )  +  H  (ti)  (Xfc  (t- )  -  %  (t-))  +  Aukjf 

(ti)  [T  (ti)  Afc  (t-)  T'^  (t,)]  “'t  (to 
[rfc  (t- )  +  H  (ti)  (SJfc  (t- )  -  %  (t-))  +  Aukj] 
q*  =  -  ^  [rfc  (t- )  +  H  (ti)  (Xfc  (t- )  -  %  (t-))  +  Arifcjj^AO^  (t” ) 

[rfc  (t“ )  +  H  (ti)  (Xfc  (t“ )  -  %  (t“))  +  Arifcj]  V  j  (3.37) 

Note  that  T  (ti)  has  been  eliminated  from  Equation  (3.37).  Since  T  (ti)  is  the  only 
variable  that  depends  upon  the  design  parameter  7,  p*^  (ti)  is  independent  of  7.  This 
means  that  Equation  (3.35)  (or,  equivalently,  Equation  (3.27))  will  yield  exactly  the  same 
hypothesis  conditional  probability  regardless  of  what  generalized  residual  vector  is  used 
in  the  GRMMAE  for  the  measurement  bias  determination  problem.  The  7  independence 
can  also  be  interpreted  as  meaning  that  /3-stripping,  described  in  Section  2.4.1,  is  always 
performed  for  applications  in  which  T  (tj)  does  not  vary  among  the  elemental  filters. 

3. 1.2  Varying  T  [ti]  In  Elemental  Filters.  As  in  the  previous  subsection,  begin  by 
substituting  the  residuals  from  the  elemental  filters  into  the  densities  in  Equation  (3.27). 
Assuming  that  each  elemental  filter  operates  on  the  same  measurement,  the  residual  charac¬ 
teristics  given  in  Equation  (3.15)  hold.  Substituting  these  characteristics  in  Equation  (3.27) 
yields 


PkiU)  = 


(27r)’"/2abs(|Tfc(ti)|)|Afc(7-)|l/2 


exp  {-^rl  {t.)Af^^  {t-  )  Vk  [t^  )  ]  pi  (ti-i) 


K 

Y _ 1 _ 

(27r)-/2abs(|T,.(7)|)|Afc(t-)|l/2 


exp 


[q*j  {ti-i) 


(3.38) 


where 
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Qj  (ti)  =  -  ]^-rf  iU)  A*  '  (ti)  r*  (U) 

=  -  \  ^k{t-)  +  U  (ti)  (xk  {t- )  -  %  (t"))  tJ  (ti)  Af'  (ti)  Tj  (ti) 

rfc(t-)  +  UHi)  (xk{t-)  - 

=  —  -  ffc  (ij  )  +  H  (tj)  (t-  )  —  Xj  (t-  Tj  (tj)  [Tj  (tj)  Afc  (tj  )  Tj  (tj)] 

Tj(ti)  rfc(t-)  +  H  (tj)  (^Sfc  (t-)  -  Sj(t“)) 

=  -l[rk{t-)  +  H(tj)(Sfc(t-)  -  Si(t“))]  A-i(t-) 

rfc  (t“ )  +  H  (tj)  (t“ )  -  Sj  (t“ ))  V  j  (3.39) 

Next,  explicitly  remove  the  component  for  the  A:*^  filter  from  the  summation  in  the 
denominator  and  factor  out  1/  (27r)™'^^  abs(|Tfc  (tj)  |)|Afc  (t“)  and  pk  (ti-i)  to  yield 


exp{gfc  (ti)} 


expjgfc  (tj)} 

V  abs(" 
j¥=k 


PjiU-i) 

Pfchi-i) 


(3.40) 


Tfc(tj) ,  Tj  (tj),  and  the  prior  probabilities  are  the  only  terms  in  Equation  (3.40) 
that  depend  upon  7.  In  order  to  maximize  p^  {ti),  7  must  be  chosen  to  minimize 


f  |Tfc  (tj)  I  \  pj  (tj_i) 

\|Tj  (ti)  I  j  Pk  (tj-i) 


(3.41) 


This  minimizes  the  contribution  of  the  incorrect  elemental  filters,  and  maximizes 
p},  (tj).  However,  setting  7/1  has  potential  to  cause  /3-dominance.  If  the  correct  filter,  k, 
is  known  a  priori  (in  which  case  there  is  no  problem  to  solve),  an  optimal  7  to  minimize 
Equation  (3.41)  can  be  found.  However,  for  a  given  problem,  this  optimal  7  may  change 
if  k  changes,  dependent  upon  the  magnitude  of  the  /3^  (tj)  and  (3*  (tj)  terms.  Figures  3.2 
and  3.3  demonstrate  this  concept.  These  figures  show  the  (3*  (tj)  terms  (including  j  =  k) 
and  mean  hypothesis  conditional  probabilities  from  1000  Monte  Carlo  simulations  of  an 
application  with  varying  Tj  (tj).  (The  specific  models  and  parameters  for  the  simulations 
are  not  presented  here  because  they  would  distract  from  the  concept  being  illustrated.)  In 
Figure  3.2,  the  correct  hypothesis  is  arbitrarily  chosen  to  be  the  hypothesis  with  the  largest 
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steady  state  P*  (ti)  value  for  the  GRMMAE  7  =  0.  For  Figure  3.2,  a  measurement  set  was 
generated  and  used  in  GRMMAFs  with  7  =  1  and  7  =  0  (i.e.,  equivalent  to  traditional  and 
post-fit  modified  MMAFs,  respectively).  From  Figure  3.2,  both  GRMMAFs  converged  to 
the  correct  hypothesis.  However,  the  GRMMAF  with  7  =  0  converged  more  quickly  than 
the  one  with  7  =  1.  Thus,  based  on  this  particular  correct  hypothesis  (i.e.,  the  elemental 
filter  with  the  largest  steady  state  P*  {tp  when  7  =  0),  the  optimal  7  value  is  0. 

Now,  consider  Figure  3.3.  The  same  application  and  system  model  was  used  for  this 
figure.  However,  the  true  parameter  value  was  changed  to  the  parameter  associated  with 
the  elemental  filter  having  the  smallest  steady  state  P*  {tp  values  when  7  =  0.  Again,  a 
measurement  set  was  generated  (associated  with  this  new  correct  parameter  value),  and 
GRMMAFs  with  7  =  1  and  7  =  0  were  used.  From  Figure  3.3,  the  GRMMAF  with 
7  =  1  (i.e.,  traditional  MMAF  equivalent)  again  converged  to  the  correct  hypothesis. 
However,  the  GRMMAE  with  7  =  0  converged  to  an  incorrect  hypothesis.  In  fact,  the 
GRMMAE  with  7  =  0  converged  to  the  same  hypothesis  as  the  GRMMAFs  in  Figure 
3.2.  This  occurred  because  the  pre-computable  Pj  {tp  term  (which  did  not  vary  between 
these  runs)  dominated  the  hypothesis  conditional  probability  calculation  in  the  GRMMAFs 
with  7  =  0.  For  the  simulations  in  Figure  3.3,  the  optimal  7  value  is  1.  Thus,  inconsistent 
optimal  7  values  are  achieved,  dependent  upon  the  relative  magnitude  of  the  beta  term 
associated  with  the  correct  hypothesis,  due  to  the  /3-dominance  effect.  The  /3-dominance 
effect  is  discussed  in  depth  in  Section  4.3.  However,  this  illustration  demonstrates  the 
/3-dominance  problem  as  it  relates  to  the  choice  of  the  optimal  7  value.  Specifically,  in 
some  cases,  the  optimal  7  artificially  induces  /3-dominance  which  causes  the  GRMMAE 
to  converge  to  one  hypothesis  regardless  of  the  realized  measurements.  Such  behavior  by 
the  GRMMAE  is  undesirable  because  it  means  that  the  GRMMAE  is  ignoring  the  real 
world  data  and  is  relying  entirely  on  its  internal  models.  Therefore,  determining  the  7 
that  maximizes  p^.  {tp  will  not  be  pursued  for  applications  where  Tj  {tp  varies  among  the 
elemental  filters. 
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Amplitude  Amplitude 


Beta  Terms  for  7=1  Average  Probability  for  7  =  1 


(a)  Beta  Terms  and  Average  Hypothesis  Conditional  Probability  for  GRMMAE  with  7  =  1 


Beta  Terms  for  7  =  0 


Average  Probability  for  7  =  0 


Measurement  Number 


Measurement  Number 


(b)  Beta  Terms  and  Average  Hypothesis  Conditional  Probability  for  GRMMAE  with  7  =  0 

Figure  3.2  Beta  Terms  and  Average  Hypothesis  Conditional  Probabilities  for  GRM- 
MAEs  with  7  =  1  and  7  =  0  and  Correct  Hypothesis  Having  Largest  Steady 
State  Beta  Term  for  7  =  0 
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Amplitude  Amplitude 


Beta  Terms  for  7=1  Average  Probability  for  7=1 


(a)  Beta  Terms  and  Average  Hypothesis  Conditional  Probability  for  GRMMAE  with  7  =  1 


Beta  Terms  for  7  =  0  Average  Probability  for  7  =  0 


(b)  Beta  Terms  and  Average  Hypothesis  Conditional  Probability  for  GRMMAE  with  7  =  0 

Figure  3.3  Beta  Terms  and  Average  Hypothesis  Conditional  Probabilities  for  GRM- 
MAEs  with  7  =  1  and  7  =  0  and  Correct  Hypothesis  Having  Smallest  Steady 
State  Beta  Term  for  7  =  0 
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3.8  Interpretation  of  f3- stripping  Using  GRMMAE  Derivation 

As  shown  in  Section  3.6,  the  GRMMAE  design  parameter,  7,  affects  only  the  fl*  {ti) 
normalization  terms.  When  /3-dominance  exists,  it  is  common  to  use  /3-stripping  to  com¬ 
pensate  for  it,  as  described  in  Section  2.4.1.  /3-stripping  is  easily  interpreted  in  the  context 
of  the  GRMMAE  derived  in  this  chapter.  Specifically,  let  7  vary  among  the  elemental  fil¬ 
ters  and  over  time.  Under  this  assumption,  each  elemental  filter  is  assigned  an  individual 
design  parameter  value  7^  {U)  at  each  measurement  epoch  ti.  Next,  choose  7j  (ti)  such  that 

abs(|T,(t,)|)  =/3,(t-)  (3.42) 

Equation  (3.31)  can  be  used  to  show 

P*{ti)  =abs(|T,(tOI)“Vi(t“)  (3.43) 

From  Equation  (3.43),  if  'Jj  (ti)  is  chosen  according  to  Equation  (3.42),  (5*  (U)  will 
equal  1,  and  the  (3*  (U)  terms  will  be  stripped  from  the  hypothesis  conditional  probability 
calculation.  Thus,  in  the  context  of  the  GRMMAE,  the  /3-stripping  method  is  interpreted 
as  choosing  jj  (ti)  such  that  Equation  (3.42)  is  satisfied.  Since,  for  linear  elemental  fil¬ 
ters,  f3*  (ti)  is  pre-computable  using  the  elemental  filter  models,  calculating  7^  (ti)  from 
Equation  (3.42)  is  also  pre-computable. 

3.9  Chapter  Summary 

This  chapter  defined  the  generalized  residual  vector.  It  gave  a  derivation  of  the  prop¬ 
erties  for  the  generalized  residual  vector  for  both  the  correct  and  incorrect  elemental  filters. 
It  described  the  effect  of  the  transformation  matrix  T  (ti).  Next,  it  provided  a  derivation 
of  the  GRMMAE’s  hypothesis  conditional  probability  formula.  Then,  the  relationship 
between  this  generalized  hypothesis  conditional  probability  (ti)  and  the  traditional  hy¬ 
pothesis  conditional  probability  pk  (ti)  was  derived.  Most  importantly,  it  demonstrated 
that  the  effect  of  the  7  design  parameter  is  to  increase  or  decrease  the  /3*  (ti)  term  mag- 
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nitude  artificially  in  the  Gaussian  densities  used  to  calculate  the  hypothesis  conditional 
probability.  Because  of  this,  the  GRMMAE  is  susceptible  to  the  /3-dominance  effect.  Fi¬ 
nally,  this  chapter  concluded  with  a  derivation  of  the  generalized  hypothesis  conditional 
probability  formula  for  applications  in  which  the  transformation  matrix  T  (tj)  varies  and 
does  not  vary  among  the  elemental  filters.  In  the  case  when  T  (tj)  does  not  vary  among 
elemental  filters,  it  demonstrated  that  the  GRMMAE  is  equivalent  to  a  traditional  MMAE, 
regardless  of  the  7  value.  In  the  case  when  T  (tj)  varies  among  the  elemental  filters,  the 
generalized  residual  hypothesis  conditional  probability  formula  remains  dependent  upon 
7  making  those  applications  susceptible  to  /3-dominance  if  7  7^  1.  Since  the  GRMMAE 
is  equivalent  to  a  traditional  MMAE  for  non- varying  T  (fj)  applications,  and  since  the 
GRMMAE  is  susceptible  to  /3-dominance  in  all  other  applications,  it  is  not  recommended 
for  most  applications.  However,  as  derived  in  this  chapter,  the  7  design  parameter  gives 
an  ability  to  exert  some  measure  of  control  over  the  magnitude  of  the  Gaussian  density 
l3j  {ti)  terms.  Therefore,  the  GRMMAE  may  be  useful  in  applications  where  /3-dominance 
is  expected.  This  is  explored  further  in  Sections  4.4  and  5.2.  Finally,  all  sections  in  this 
chapter  are  research  contributions  of  this  dissertation. 
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IV.  Simulations  and  Analysis 


4-1  Residual  Properties  Verifieation 

Sections  4.1.2  through  4.1.4  present  examples  to  verify  the  generalized  residual  vector 
and  generalized  residual  covariance  matrix  equations  derived  in  Chapter  III.  These  exam¬ 
ples  are  presented  in  the  same  order  as  the  derivations  in  Sections  3.3.1  and  3.3.4  with 
which  they  correspond.  The  examples  all  use  a  common  truth  model  which  is  described 
in  Section  4.1.1.  Since  GPS  carrier- phase  integer  ambiguity  resolution  is  the  application 
that  inspired  the  GRMMAE  development,  the  truth  model  is  designed  to  mimic  that 
application.  However,  in  an  effort  to  keep  the  GRMMAE  demonstration  general,  GPS 
carrier- phase  integer  ambiguity  resolution  is  not  specifically  used  in  these  sections.  Ambi¬ 
guity  resolution  will  be  used  in  Section  4.2  to  demonstrate  another  GRMMAE  property.  In 
each  of  these  examples,  a  different  parameter  is  assumed  to  be  unknown,  and  a  GRMMAE 
is  built  to  determine  the  correct  parameter  value  given  in  the  common  truth  model. 

4.1.1  Truth  Model.  The  truth  model  for  the  examples  in  Sections  4.1.2  through 
4.1.4  is  a  simple  1-D  system.  It  describes  a  vehicle  moving  along  a  line.  The  system 
state  vector  consists  of  position,  velocity,  and  acceleration  states,  all  along  the  line.  The 
dynamics  model  is  a  first  order  Gauss-Markov  acceleration  (EOGMA)  model  given  by 


x(t)  =  F(t)x(t)  -F  w(t) 


(4.1) 


where  x  (t)  is  the  state  vector  consisting  of  the  position,  velocity,  and  acceleration  states 


x(t) 


P(t) 

v(t) 

a(t) 


and  F  (t)  is  the  dynamics  matrix  given  by 


(4.2) 
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0  1 


0 


(4.3) 


F(t)  =  0  0  1 

_  0  0  -1/ta  _ 

where  ta  is  the  acceleration  time  constant.  Finally,  in  Equation  (4.1),  w  (t)  is  the  contin¬ 
uous  time  dynamics  noise  given  by 

0 

w(t)  =  0  (4.4) 

Wa  (t) 

where  Wa  (t)  is  the  zero-mean,  white,  Gaussian  noise  driving  the  acceleration  state.  The 
covariance  kernel  of  the  noise  vector  is  given  by 

F[w(t)w^(r)]  =  Q  (t)  (5  (t  —  r)  (4.5) 

The  FOGMA  model  is  commonly  used  in  navigation  applications  including  GPS 
carrier- phase  integer  ambiguity  resolution.  It  assumes  that  the  vehicle  acceleration  is  well 
described  with  a  first-order  lag  driven  by  white,  Gaussian  noise.  From  Equations  (4.1) 
through  (4.4),  the  model  for  the  acceleration  state  is  given  by 

^  a(t)  =  a(t)  +  Wa  (t) 

where  a  (t)  is  the  vehicle  acceleration  in  1-D  and  Wg  (t)  is  the  continuous-time  white,  Gaus¬ 
sian  noise  driving  the  model.  Ghanging  the  acceleration  time  constant  and/or  the  strength, 
Q  (t),  of  the  driving  noise  changes  the  vehicle  motion  described  by  the  truth  model.  A  true 
value  of  ta  =  2  is  arbitrarily  chosen  for  these  examples.  The  strength  of  Wa  (t)  is  used  in 
the  calculation  of  the  equivalent  discrete-time  noise  covariance  matrix  Qrf(ti),  which  will 
be  discussed  shortly. 
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Since  the  dynamics  matrix  in  Equation  (4.3)  is  time-invariant,  the  state  transition 
matrix  is  calculated  using  the  matrix  exponential  [40] 

(4.6) 

where  At  =  ti  —  U-i. 

Using  ta  =  2  and  Equations  (4.3)  and  (4.6),  the  true  state  dynamics  model  is 

1  1  0.426 

x(U+i)  =  0  1  0.787  x(U)  Wrf(U)  (4.7) 

0  0  0.607 

In  Equation  (4.7),  (ti)  is  the  equivalent  discrete-time  dynamics  driving  noise  calculated 
from  the  continuous-time  noise  w(t).  The  equivalent  discrete-time  noise  is  zero-mean, 
white,  Gaussian  with  covariance  of  (ti)  calculated  from  Q  (t)  [29] 

ti+i 

Qd(ti)  =  j  $(U+i,r)Q(r)$'^(U+i,r)(ir  (4.8) 

U 

Two  independent  measurement  sensors  are  used  for  these  examples.  Each  sensor 
measures  the  1-D  position  of  the  vehicle  along  a  line.  The  measurement  model  is  given  by 

1  0  0 

z(ti)  =  x(ti)  +  \j(ti)  (4.9) 

1  0  0 

where  v  (ti)  is  the  measurement  noise  vector. 

The  equivalent  discrete-time  dynamics  driving  noise  and  measurement  noise  are  in¬ 
dependent,  zero-mean,  white,  Gaussian  noises  with  covariances  of  (U)  and  R(U),  re¬ 
spectively.  Two  noise  scenarios  are  used,  a  High  Q^/Low  R  scenario  (which  is  similar  to 
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the  GPS  carrier-phase  ambiguity  resolution  application)  and  a  Low  Q^/High  R  scenario. 
These  relative  values  are  chosen  to  create  scenarios  in  which  the  elemental  filters  emphasize 
one  model  (dynamics  or  measurement)  over  the  other.  The  noise  covariance  matrices  for 
these  scenarios  are 


High  Qf^/Low  R  Scenario 


Low  Q^/High  R  Scenario 


Qd  {ti) 


R(L) 


0.0454 

0.0838 

0.0646 

0.0454 

0.0838 

0.0646 

0.0838 

0.1548 

0.1193 

Qd(ti)  = 

0.0838 

0.1548 

0.1193 

0.0646 

0.1193 

0.0920 

0.0646 

0.1193 

0.0920 

0.0025 

0 

R(L)  = 

25 

0 

0 

0.0025 

0 

25 

(4.10) 


The  initial  conditions  for  each  scenario  are 


0 

25 

0 

0 

1 

Po  = 

0 

100 

0 

0 

0 

0 

10 

(4.11) 


Finally,  in  each  example  that  follows  (including  the  measurement  bias  example)  the 
true  measurements  are  unbiased.  Thus,  in  the  measurement  bias  example,  the  true  bias  is 
given  by 


^true 


0 

0 


(4.12) 


The  truth  model  given  in  Equations  (4.1)  through  (4.12)  is  used  for  the  examples 
in  the  next  four  sections.  These  examples  demonstrate  the  generalized  residual  mean  and 
covariance  equations  derived  in  Sections  3.3.1  through  3.3.4.  The  examples  are  presented 
in  the  same  order  as  the  derivations  in  Chapter  III.  The  GRMMAEs  applied  to  this  truth 
model  vary  according  to  the  unknown  parameter  being  estimated.  The  GRMMAE  models 
for  each  example  are  given  in  the  corresponding  section. 
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4-1-2  Example:  Unknown  Measurement  Bias  Seenario-  This  example  is  used  to 
demonstrate  the  residual  properties  derived  in  Section  3.3.1.  The  truth  model  presented 
in  Section  4.1.1  is  used  for  this  example.  The  elemental  filters  in  the  GRMMAE  are 
identical  but  act  upon  corrected  measurements  given  by  Equation  (3.7),  repeated  here  for 
convenience 


(^i)  —  ^  (li) 


(3.7) 


where  z  (tj)  is  the  biased  sensor  measurement,  and  is  the  hypothesized  bias  for  the 
elemental  filter.  Thus,  the  elemental  filter  model  for  this  example  is  given  by  Equations 
4.1  through  4.12.  Three  elemental  filters  are  used,  hypothesizing  different  measurement 
biases,  which  is  very  similar  to  the  carrier-phase  integer  ambiguity  resolution  approach  in 
[17,  18].  The  hypothesized  measurement  biases  are 


0 

5 

0 

m  = 

ns  = 

0 

10 

5 

(4.13) 


where  nj  is  the  hypothesized  bias  vector  implemented  in  the  elemental  filter  (j  =  1,2,3). 
The  true  measurement  bias  is  ni. 

This  system  was  simulated  1000  times,  running  each  simulation  for  100  measurement 
epochs.  Sample  averages  of  the  generalized  residual  vectors  and  the  generalized  residual 
covariance  matrices  were  calculated  at  each  measurement  epoch.  Thus,  the  1000  Monte 
Carlo  runs  were  used  to  yield  a  100-step  time  history  of  2-dimensional  mean  vectors  (2 
measurement  sources,  100  measurement  epochs)  and  100  2  x  2  sample  covariance  matrices 
(one  for  each  epoch).  Additionally,  sample  averages  of  the  state  vectors  prior  to  mea¬ 
surement  incorporation  were  calculated  for  each  elemental  filter.  Thus,  100  state  vectors 
(one  for  each  epoch),  Xj  (t~),  were  produced  for  each  of  the  3  elemental  filters.  Since  the 
H  (tj)  and  K  (tj)  matrices  do  not  vary  from  one  Monte  Carlo  run  to  the  next,  it  is  not 
necessary  to  calculate  sample  averages  of  H  (tj)  ,  K  (t*) ,  or  T  (tj).  Additionally,  since  the 
hypothesized  measurement  biases  do  not  change  from  one  Monte  Carlo  run  to  the  next,  it 
is  not  necessary  to  compute  sample  averages  of  the  difference  in  the  hypothesized  biases 
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Arifcj.  This  process  was  repeated  a  total  of  4  times,  once  for  each  combination  of  the  noise 
covariance  matrices  in  Equation  (4.10)  with  7  equal  to  both  0  and  1.  This  is  depicted 
graphically  by  the  “X”  marks  in  Table  4.1 

For  each  Monte  Carlo  simulation  listed  in  Table  4.1,  the  sample  averages  were  used 
to  calculate  the  mean  of  the  generalized  residual  vector  for  the  incorrect  elemental  filters  (2 
and  3)  using  Equation  (3.11),  a  modified  portion  of  which  is  repeated  here  for  convenience 

uir*  (ti)  =  T  (ti)  [H  (ti)  {2k,,,,  {t~ )  -  {f-))  +  Arifcj]  (3.11) 

J theory 

where  2k,,„  (t~)  ,  and  2j„„  are  the  sample  averages  taken  from  the  1000  Monte 
Carlo  runs,  and  the  subscripts  aver  and  theory  denote  sample  averages  and  a  residual 
mean  vector  calculated  using  the  theory  from  Chapter  III,  respectively.  It  was  necessary 
to  calculate  sample  averages  of  the  state  vector  estimates,  2k„„  (t~)  and 
because  the  true  state  varied  in  each  Monte  Carlo  run  due  to  the  dynamics  model  noise. 
The  subscripts  aver  and  theory  were  added  here  to  distinguish  between  realized  random 
variable  values  in  the  1000  Monte  Carlo  simulations  {aver)  and  values  calculated  using  the 
equations  in  Chapter  III  (theory). 

Figure  4.1  shows  a  plot  of  the  theoretical  generalized  residuals,  calculated  using  the 
modified  version  of  Equation  (3.11),  and  the  realized  sample  average  residual  vectors  for 
elemental  filters  2  and  3  (the  incorrect  filters).  In  this  figure,  the  Low  Q^/High  R  noise 
scenario  was  used.  For  that  scenario,  the  elemental  Kalman  filters  place  increased  emphasis 
on  the  filter  dynamics  model  and  less  emphasis  on  the  measurements.  Additionally,  the 
traditional  residual  (7  =  1)  was  used.  Finally,  in  order  to  make  the  figure  more  clear,  only 

Table  4.1  Matrix  of  Monte  Carlo  Test  Scenarios:  Combinations  of  Noise  Matrices  and 
7  Values 


Generalized  Re 
Traditional  (7  =  1) 

isidual  Type 
Post-Fit  (7  =  0) 

Noise  High  Q^^/Low  R 

X 

X 

Scenario  Low  Q^j/High  R 

X 

X 
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Theoretical  and  Sample  Average  Generalized  Residuals 


(a)  Elemental  Filter  2 


Theoretical  and  Sample  Average  Generalized  Residuals 


(b)  Elemental  Filter  3 


Figure  4.1  Theoretical  and  Sample  Average  Generalized  Residuals  from  1000  Simula¬ 
tions  for  Measurement  Source  1  in  Low  (^O/High  R  Noise  Scenario  with 
Traditional  Residual  (7  =  1)-  Measurement  Bias  Example 
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the  residuals  associated  with  measurement  source  1  are  shown.  The  plot  for  measurement 
source  2  was  similar.  Figure  4.1  demonstrates  that  the  sample  mean  generalized  residuals 
track  the  theoretical  mean  generalized  residuals  from  Equation  (3.11).  Relative  to  the 
generalized  residual  covariance,  the  realized  sample  mean  closely  matches  the  theoretical 
mean.  The  simulations  for  the  other  combinations  of  noise  covariance  matrices  and  gamma 
values  in  Table  4.1  produced  results  similar  to  Figure  4.1  and  are  not  shown.  Figure  4.1 
demonstrates  that  Equation  (3.11)  correctly  calculated  the  mean  value  for  the  generalized 
residuals  from  the  incorrect  elemental  filters,  based  on  the  traditional  residual  from  the 
correct  elemental  filter. 

Next,  the  covariance  of  the  generalized  residuals  from  the  incorrect  elemental  fil¬ 
ters  was  confirmed.  For  the  measurement  bias  application,  this  formula  is  given  by 
Equation  (3.13),  a  modified  version  of  which  is  repeated  here  for  convenience 

=T{t.)Ai,„,((-)T^(*i)  (3.13) 

where  the  subscript  theory  indicates  a  matrix  calculated  using  the  theory  in  Chapter  III. 
The  new  subscript  filt  indicates  that  the  elemental  filter-computed  covariance  was  used. 
The  filter-computed  residual  covariance  was  given  immediately  after  Equation  (2.47)  and 
is 


(f-)  =  H  iU)  Pk  {t-)  iU)  +  R  iU)  (4.14) 

None  of  the  matrices  above  vary  from  one  Monte  Carlo  run  to  the  next.  Therefore,  there 
is  no  need  to  calculate  a  sample  average  value  for  the  elemental  filter-computed  covariance 
matrix. 

Figure  4.2  shows  a  representative  plot  of  the  realized  sample  covariance  matrix  and 
the  theoretical  covariance  matrix  calculated  from  Equation  (3.13).  Since  there  were  two 
measurement  sources,  the  residual  covariance  matrices  were  2x2.  In  Figure  4.2,  the 
line  marked  “Theoretical  Variance  for  Sensors  1  and  2”  is  a  plot  of  the  variance  values 
from  the  diagonal  of  the  covariance  matrix.  In  this  case,  the  two  variance  values  are  the 
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Theoretical  and  Sample  Generalized  Residual  Covariance  Matrix  Elements 


(a)  Elemental  Filter  2 


Theoretical  and  Sample  Generalized  Residual  Covariance  Matrix  Elements 


(b)  Elemental  Filter  3 


Figure  4.2  Theoretical  and  Sample  Average  Generalized  Residual  Covariance  Matrix 
Elements  from  1000  Simulations  in  Low  Q(i/High  R  Noise  Scenario  with 
Post-Fit  Residuals  (7  =  0)-  Measurement  Bias  Example 
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same  and  only  one  line  appears  in  the  plot,  representing  both  sensors.  The  line  marked 
“Theoretical  Covariance  for  Sensors  1  and  2”  is  a  plot  of  the  off-diagonal  elements  of  the 
covariance  matrix.  Since  covariance  matrices  are  symmetric,  there  is  only  a  need  to  plot 
a  single  line  for  the  covariance  values.  The  plots  of  “Realized  Variance”  and  “Realized 
Covariance”  are  the  values  from  the  realized  sample  covariance  matrices  corresponding 
to  the  theoretical  values.  Figure  4.2  shows  the  variance/covariance  values  for  the  Low 
Qrf/High  R  noise  scenario  using  post-fit  residuals  (7  =  0).  Figure  4.2  demonstrates  that 
the  realized  sample  variance  and  covariance  values  track  the  corresponding  theoretical 
values  with  errors  that  are  an  order  of  magnitude  smaller  than  the  theoretical  values.  This 
behavior  was  true  for  all  combinations  of  noise  covariance  matrices  and  gamma  values  in 
Table  4.1,  and  indicates  that  Equation  (3.13)  correctly  calculated  the  generalized  residual 
covariance  matrix.  However,  Figure  4.2  is  shown  because  it  demonstrates  unexpected 
behavior.  Typically,  the  residual  variance  is  large  at  first  and  settles  down  to  a  smaller 
steady  state  value.  In  Figure  4.2,  though,  the  residual  variance  values  start  smaller  and 
increase  to  their  steady  state  value.  This  was  caused  by  the  transformation  matrix  T  (tj). 
This  unexpected  behavior  will  be  discussed  in  Chapter  V  in  relation  to  a  recommendation 
for  future  research  into  use  of  the  GRMMAE  for  /3-dominance  compensation. 

Einally,  based  on  Equations  (3.10)  and  (3.13),  for  this  application,  the  generalized 
residuals  from  the  incorrect  elemental  filters  equal  the  generalized  residuals  from  the  correct 
elemental  filter,  offset  by  a  bias.  Eigure  4.3  demonstrates  this  property.  In  this  figure,  the 
generalized  residual  associated  with  measurement  source  1  from  all  three  elemental  filters  is 
shown.  Unlike  previous  figures,  a  single  Monte  Carlo  run  was  used  for  this  demonstration. 
This  ensures  that  a  direct  comparison  can  be  made  between  the  elemental  filters.  If  a 
sample  average  had  been  used,  it  is  possible  (though  not  likely)  to  have  multiple  runs  which 
individually  exhibit  incorrect  behavior  while  still  maintaining  correct  average  behavior.  By 
using  a  single  run,  this  possibility  was  avoided.  Eigure  4.3  shows  the  generalized  residuals 
for  the  High  Q^^/Low  R  noise  covariance  scenario  with  post-fit  residuals  (i.e.,  7  =  0).  In 
Eigure  4.3,  the  generalized  residuals  from  elemental  filters  2  and  3  overlap  due  to  the  choice 
of  measurement  bias  hypotheses,  and  thus,  only  one  line  appears  in  the  plot  representing 
both  filters.  A  comparison  of  the  residuals  from  elemental  filter  1  (the  correct  hypothesis) 


4-10 


Generalized  Residuals  for  Measurement  Sensor  1 


Figure  4.3  Single  Run  Generalized  Residual  for  Measurement  Sensor  1  in  High  Q^^/Low 
R  Noise  Scenario  with  Post-Fit  Residuals  (7  =  0) 

with  those  from  elemental  filters  2  and  3  clearly  shows  that  the  residuals  from  the  incorrect 
elemental  filters  equal  the  residuals  from  the  correct  elemental  filter,  with  an  added  bias. 

4.1.3  Example:  Unknown  Scenario.  This  example  demonstrates  an  applica¬ 
tion  with  an  unknown  equivalent  discrete-time  dynamics  noise  covariance  matrix.  Equa¬ 
tions  to  calculate  the  generalized  residuals  for  the  incorrect  elemental  filters  from  the  tra¬ 
ditional  residuals  for  the  correct  filter  are  given  in  Section  3.3.2.  As  with  other  examples, 
the  truth  model  presented  in  Section  4.1.1  was  used  for  this  example. 

An  example  such  as  this  might  arise,  for  example,  in  a  target  tracking  problem.  In 
this  application,  the  FOGMA  model  in  Equations  (4.1)  through  (4.7)  adequately  describes 
the  dynamics  of  a  vehicle.  However,  changes  in  target  dynamics  are  described  by  varying 
the  strength  of  the  continuous-time  dynamics  noise  w{t).  Eor  instance,  a  target  at  rest 
would  be  represented  with  a  very  small  dynamics  model  noise  strength.  Various  dynamic 
operating  modes  (accelerating,  constant  velocity,  etc.)  would  be  represented  with  other 
dynamics  noise  strengths.  Those  varying  values  of  strength  for  w  (f)  lead  to  varying  {ti) 
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for  each  operating  mode.  The  GRMMAE  is  then  used  to  determine  which  hypothesized 
Qd  (ti)  best  matches  the  motion  of  the  vehicle,  allowing  vehicle  tracking. 

In  this  example,  three  elemental  filters  were  used,  hypothesizing  three  values  for 
the  equivalent  discrete-time  noise  covariance.  The  hypothesized  values  for  (ti)  were 
arbitrarily  chosen  to  be 


Qditi)  2.5Qd{ti)  6.7Qd{U) 

where  (ti)  is  the  value  given  in  Equation  (4.10).  Two  noise  scenarios  and  two  values  of 
7  were  used  as  shown  in  Table  4.1,  for  a  total  of  four  sets  of  Monte  Carlo  simulations.  Each 
scenario  defined  in  Table  4.1  was  simulated  1000  times  with  100  measurement  epochs  used 
in  each  Monte  Carlo  simulation.  As  in  the  measurement  bias  example,  sample  averages  of 
the  variables  necessary  to  confirm  the  Equations  in  Section  3.3.2  were  calculated.  Since 
the  H  (tj)  and  Kj  (tj)  matrices  do  not  vary  from  one  Monte  Carlo  run  to  the  next,  it  was 
not  necessary  to  calculate  sample  averages  of  H  (ti)  ,  Kj  (U)  ,  or  Tj  (ti).  Notice,  however, 
that  unlike  the  previous  example,  Kj  (tj)  and  Tj  (tj)  do  vary  among  the  elemental  filters. 

Eor  each  Monte  Carlo  simulation  listed  in  Table  4.1,  the  sample  averages  were  used 
to  calculate  the  mean  of  the  generalized  residual  vector  for  the  incorrect  elemental  filters  (2 
and  3)  using  Equation  (3.15),  a  modified  portion  of  which  is  repeated  here  for  convenience 

mr*  {ti)  =  Tj  (ti)  H  (ti)  {f- )  -  {f- )]  (3.15) 

^theory 

where  ^kaver  (^i~)  ’  ^javer  (^7 )  the  Sample  averages  taken  from  the  1000  Monte 

Carlo  runs,  and  the  subscripts  aver  and  theory  denote  sample  averages  and  a  residual 
mean  vector  calculated  from  the  theory  in  Chapter  III,  respectively. 

Eigure  4.4  shows  a  plot  of  the  theoretical  generalized  residuals,  calculated  using  the 
modified  version  of  Equation  (3.15),  and  the  realized  sample  average  residual  vectors  for 
elemental  filters  2  and  3  (the  incorrect  filters).  In  this  figure,  the  Low  Q^/High  R  noise 
scenario  was  used.  Eor  that  scenario,  the  elemental  Kalman  filters  place  increased  empha¬ 
sis  on  the  filter  dynamics  model  and  less  emphasis  on  the  measurements.  Additionally, 


4-12 


Measurement  Number 


(a)  Elemental  Filter  2 


Theoretical  and  Sample  Average  Generalized  Residuals 


(b)  Elemental  Filter  3 


Figure  4.4  Theoretical  and  Sample  Average  Generalized  Residuals  from  1000  Simula¬ 
tions  for  Measurement  Source  1  in  Low  Qrf/High  R  Noise  Scenario  with 
Post-Fit  Residual  (7  =  0)-  Unknown  Example 
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the  post-fit  residual  (7  =  0)  was  used.  Finally,  in  order  to  make  the  figure  more  clear, 
only  the  residuals  associated  with  measurement  source  1  are  shown.  Figure  4.4  shows  that, 
relative  to  the  generalized  residual  standard  deviation,  the  realized  sample  mean  gener¬ 
alized  residual  exhibited  a  good  match  to  the  theoretical  mean.  This  demonstrates  that 
Equation  (3.15)  correctly  calculated  the  mean  generalized  residual  vector  for  the  incorrect 
elemental  filters  based  on  the  traditional  residual  vector  from  the  correct  elemental  filter. 

Next,  the  covariance  of  the  generalized  residuals  from  the  incorrect  elemental  fil¬ 
ters  was  confirmed.  For  the  unknown  application,  this  formula  is  also  given  by 
Equation  (3.15),  a  modified  version  of  which  is  repeated  here  for  convenience 

77.,,  (*.)  =  7  (ti)  ((-)  Tj  {(,)  (3.15) 

where  the  subscripts  theory  and  filt  indicate  a  covariance  matrix  calculated  from  the 
theory  in  Chapter  III  and  the  elemental  filter-computed  covariance  matrix,  respectively. 
As  before,  the  filter-computed  covariance  matrix  does  not  change  from  one  run  to  the  next. 
Therefore,  there  is  no  need  to  calculate  a  sample  average  value  for  it. 

Figure  4.5  shows  a  representative  plot  of  the  realized  sample  covariance  matrix  and 
the  covariance  matrix  calculated  from  Equation  (3.15).  In  Figure  4.5,  the  line  marked 
“Theoretical  Variance  for  Sensors  1  and  2”  is  a  plot  of  the  variance  values  from  the  diagonal 
of  the  covariance  matrix.  In  this  case,  the  two  variance  values  were  the  same  and  only  1  line 
appears  in  the  plot,  representing  both  sensors.  The  line  marked  “Theoretical  Covariance 
for  Sensors  1  and  2”  is  a  plot  of  the  off-diagonal  elements  of  the  covariance  matrix.  The 
plots  of  “Realized  Variance”  and  “Realized  Covariance”  are  the  values  from  the  realized 
sample  covariance  matrices  corresponding  to  the  theoretical  values.  Figure  4.5  shows 
the  variance /covariance  values  for  the  Low  Q^^/High  R  noise  scenario  using  traditional 
residuals  (7  =  1).  Figure  4.5  demonstrates  that  the  realized  sample  variance  and  covariance 
values  track  the  corresponding  theoretical  values  with  errors  that  are  an  order  of  magnitude 
smaller  than  the  theoretical  values.  This  behavior  was  true  for  all  combinations  of  noise 
covariance  matrices  and  gamma  values  in  Table  4.1,  and  indicates  that  Equation  (3.15) 
correctly  calculated  the  generalized  residual  covariance  matrix. 
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Theoretical  and  Sample  Generalized  Residual  Covariance  Matrix  Elements 


(a)  Elemental  Filter  2 


Theoretical  and  Sample  Generalized  Residual  Covariance  Matrix  Elements 


(b)  Elemental  Filter  3 


Figure  4.5  Theoretical  and  Sample  Generalized  Residual  Covariance  Matrix  Elements 
from  1000  Simulations  in  Low  Qrf/High  R  Noise  Scenario  with  Traditional 
Residuals  (7  =  1)-  Unknown  Example 
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4-1-4  Example:  Unknown  F  (t)  Scenario.  This  example  confirms  the  equations 
for  the  generalized  residual  vector  and  covariance  matrix,  derived  in  Section  3.3.3,  for 
an  application  in  which  the  matrix  F  (t)  describing  the  homogeneous  system  dynamics 
is  subject  to  uncertainty.  This  application  could  arise  in  a  target  identification  problem. 
In  that  application,  various  targets  would  be  described  by  the  general  FOGMA  model 
given  in  Equation  (4.3).  Each  hypothesized  target  would  have  a  different  acceleration  time 
constant  ta,  which  in  turn  leads  to  each  hypothesis  having  a  different  F  (t)  matrix.  A 
target  being  tracked  would  then  be  identified  by  matching  the  observed  motion  to  the  best 
hypothesis  for  F  (t) . 

In  this  example,  three  elemental  filters  were  used,  hypothesizing  three  values  for  the 
acceleration  time  constant  ta-  The  hypothesized  values  for  ta  were  arbitrarily  chosen  to  be 
2,  4,  and  8.  Substitution  of  these  acceleration  constants  into  Equation  (4.3),  yields  the 
three  hypotheses  for  F  (t) 


0 

1 

0 

0 

1 

0 

0 

1 

0 

0 

0 

1 

F2(t)  = 

0 

0 

1 

F3(t)  = 

0 

0 

1 

0 

0 

-0.5 

0 

0 

-0.25 

0 

0 

-0.125 

(4.15) 


From  Section  4.1.1,  Fi  (t)  was  the  true  value  of  the  homogeneous  dynamics  matrix. 

The  analysis  steps  for  the  1000  Monte  Carlo  simulations  in  Section  4.1.3  were  re¬ 
peated  for  this  application,  with  similar  results  and  conclusions.  Figures  4.6  and  4.7 
demonstrate  that  Equation  (3.15)  correctly  calculated  the  mean  vectors  and  covariance 
matrices  of  the  generalized  residuals  from  the  incorrect  elemental  filters.  Each  figure  is 
labelled  according  to  the  noise  covariance  scenario  and  type  of  generalized  residual  used. 

4.1.5  Example:  Unknown  H  Scenario.  This  example  demonstrates  an  application 
with  an  unknown  measurement  noise  covariance  matrix.  The  equations  for  the  generalized 
residual  and  its  covariance  for  such  an  application  are  discussed  in  Section  3.3.4.  This 
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Figure  4.6 


Theoretical  and  Sample  Average  Generalized  Residuals 
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Theoretical  and  Sample  Average  Generalized  Residuals  from  1000  Simula¬ 
tions  for  Measurement  Source  1  in  High  Q^^/Low  R  Noise  Scenario  with 
Post-Fit  Residual  (7  =  0)-  Unknown  F  (t)  Example 
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Theoretical  and  Sample  Generalized  Residual  Covariance  Matrix  Elements 
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Figure  4.7  Theoretical  and  Sample  Generalized  Residual  Covariance  Matrix  Elements 
from  1000  Simulations  in  High  Q^/Low  R  Noise  Scenario  with  Post-Fit  Resid¬ 
uals  (7  =  0)-  Uknown  F  (t)  Example 
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section  presents  results  to  confirm  the  validity  of  those  equations.  As  with  other  examples, 
the  truth  model  presented  in  Section  4.1.1  was  used  for  this  example. 

An  example  of  such  an  application  from  the  navigation  field  is  navigation  using  a 
GPS  aided  Inertial  Navigation  System  (INS)  in  a  GPS  jamming  environment.  As  jamming 
strength  increases,  the  accuracy  of  the  GPS  receiver  decreases.  This  is  modelled  by  in¬ 
creasing  the  measurement  noise  covariance  R  (tj).  If  the  jamming  becomes  strong  enough, 
GPS  signal  tracking  is  lost  altogether,  and  the  INS  must  run  without  GPS  updates,  which 
is  modelled  by  another  measurement  noise  covariance  level.  A  GRMMAE  that  can  de¬ 
termine  the  appropriate  measurement  noise  covariance  matrix  would  allow  the  navigation 
system  to  continue  operation  through  a  degradation  and  eventual  loss  of  the  GPS  sensor. 
Additionally,  it  would  automatically  account  for  the  re-introduction  of  GPS  as  the  system 
left  the  jamming  field. 

In  this  example,  three  elemental  filters  were  used,  hypothesizing  three  values  for  the 
measurement  noise  covariance  matrix.  The  hypothesized  values  for  R  were  arbitrarily 
chosen  to  be 

R(ti)  2.5R(ti)  6.7R(ti) 

where  R(tj)  is  the  value  given  in  Equation  (4.10). 

The  analysis  steps  for  the  1000  Monte  Garlo  simulations  in  Section  4.1.3  were  re¬ 
peated  for  this  application,  with  similar  results  and  conclusions.  Eigures  4.8  and  4.9 
demonstrate  that  Equation  (3.15)  correctly  calculated  the  mean  vectors  and  covariance 
matrices  of  the  generalized  residuals  from  the  incorrect  elemental  filters.  Each  figure  is 
labelled  according  to  the  noise  covariance  scenario  and  type  of  generalized  residual  used. 

4-1.6  Measurement  Bias  Parameter  with  7  =  0.5.  In  the  previous  four  examples, 
the  GRMMAE  design  parameter  7  was  always  set  to  0  or  1  (post-fit  or  traditional  residuals, 
respectively).  This  was  done  arbitrarily  because  these  two  values  make  the  most  physical 
sense.  However,  at  no  point  in  the  GRMMAE  derivation  was  7  constrained  to  be  0  or  1. 
In  fact,  7  is  allowed  to  be  any  real,  scalar  value.  Therefore,  the  measurement  bias  example 
was  repeated  for  the  High  Q^/Low  R  noise  scenario  of  Equation  (4.10)  with  7  =  0.5 
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Theoretical  and  Sample  Average  Generalized  Residuals 


(a)  Elemental  Filter  2 


Theoretical  and  Sample  Average  Generalized  Residuals 


(b)  Elemental  Filter  3 

Figure  4.8  Theoretical  and  Sample  Average  Generalized  Residuals  from  1000  Simula¬ 
tions  for  Measurement  Source  1  in  High  Q^^/Low  R  Noise  Scenario  with 
Traditional  Residual  (7  =  1)-  Unknown  R  Example 
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Theoretical  and  Sample  Generalized  Residual  Covariance  Matrix  Elements 
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Theoretical  and  Sample  Generalized  Residual  Covariance  Matrix  Elements 
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Figure  4.9  Theoretical  and  Sample  Generalized  Residual  Covariance  Matrix  Elements 
from  1000  Simulations  in  Low  Q^/High  R  Noise  Scenario  with  Traditional 
Residuals  (7  =  1)-  Unknown  R  Example 
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(i.e.,  a  linear  combination  of  the  post-fit  and  traditional  residuals).  The  results  of  these 
1000  Monte  Carlo  runs,  shown  in  Figures  4.10  and  4.11,  were  similar  to  those  in  previous 
sections,  demonstrating  that  Equations  (3.11)  and  (3.13)  for  the  mean  and  covariance  of 
r*  (ti)  hold  for  7  =  0.5,  as  expected. 

4-2  GRMMAE  for  GPS  Garrier-Phase  Integer  Ambiguity  Resolution 

Carrier-phase  integer  ambiguity  resolution  is  an  example  of  a  GRMMAE  application 
with  a  constant  measurement  bias  parameter.  This  general  class  of  problems  is  presented 
in  Sections  3.3.1  and  3.7.1.  Previous  research  demonstrated  that  a  GRMMAE  with  7  =  0 
(i.e.,  using  post-fit  residuals)  will  successfully  resolve  GPS  carrier-phase  integer  ambiguities 
[17,  18].  A  research  contribution  of  this  dissertation  is  the  derivation  in  Section  3.7.1 
which  shows  that  the  hypothesis  conditional  probabilities  calculated  for  this  application 
are  independent  of  the  choice  of  7.  Thus,  while  the  use  of  post-fit  residuals  is  valid  for 
this  application,  the  use  of  traditional  residuals  is  equally  valid.  This  section  will  use  GPS 
carrier-phase  measurement  data  to  demonstrate  that  GRMMAE’s  with  7  =  0  and  7  =  1 
produce  equivalent  results  when  the  parameter  of  interest  is  an  unknown  measurement 
bias. 

Before  proceeding,  it  is  important  to  recognize  that  the  GRMMAE  for  GPS  carrier- 
phase  integer  ambiguity  resolution  does  not  meet  the  linearity  assumption  in  Sections  3.3.1 
and  3.7.1.  In  those  sections,  it  is  assumed  that  the  GRMMAE  elemental  filters  are  linear 
Kalman  filters.  However,  the  measurement  model  for  GPS  is  a  nonlinear  model,  leading  to 
the  use  of  extended  Kalman  filters  in  the  GRMMAE  (recall  that  use  of  extended  Kalman 
filters  is  a  non-optimal,  ad  hoc  extension  to  the  MMAE  and  GRMMAE).  Therefore,  the 
measurement  matrices  Hj  (ti)  in  the  elemental  filters  are  first-order  approximations  as 
given  in  Equation  (2.26)  and  Table  2.2.  Because  the  measurement  matrices  are  linearized 
around  the  current  state  estimate,  and  each  elemental  filter  produces  a  unique  position 
estimate,  the  measurement  matrices  in  the  elemental  filters  will  not  be  equal  in  all  elemental 
filters  as  assumed  in  Section  3.7.1.  This,  in  turn,  leads  directly  to  the  transformation 
matrices  Tj  (tj)  varying  in  the  elemental  filters,  which  violates  the  key  assumption  that 
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Figure  4.10 
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Theoretical  and  Sample  Average  Generalized  Residuals  from  1000  Simula¬ 
tions  for  Measurement  Source  1  in  High  Q^/Low  R  Noise  Scenario  with 
7  =  0.5  -  Measurement  Bias  Example 
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Figure  4.11 


Theoretical  and  Sample  Generalized  Residual  Covariance  Matrix  Elements 
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Theoretical  and  Sample  Generalized  Residual  Covariance  Matrix  Elements 
from  1000  Simulations  in  High  Q^/Low  R  Noise  Scenario  with  7  =  0.5  - 
Measurement  Bias  Example 
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Tj  (ti)  =  Tfc  (ti)  for  all  j  and  k.  However,  because  of  measurement  geometry,  this  effect 
is  not  significant  for  the  GPS  carrier-phase  integer  ambiguity  resolution  application. 

When  performing  GPS  carrier-phase  integer  ambiguity  resolution,  the  measurement 
matrices  are  composed  of  unit  line  of  sight  vectors  pointing  from  the  estimated  user  po¬ 
sition  toward  the  satellites  [39].  However,  the  user  position  estimates  (produced  using 
the  hypothesized  ambiguity  vectors)  will  be  separated  by  short  distances  on  the  order  of 
centimeters.  Since  the  GPS  satellites  are  tens  of  millions  of  meters  away  from  the  user 
and  the  hypothesized  user  positions  are  separated  by  centimeters,  the  differences  in  the 
unit  line  of  site  vectors  (and,  hence  the  measurement  matrices)  are  so  small  that  they  can 
be  neglected.  Thus,  while  the  measurement  matrices  are  not  mathematically  equal  in  all 
elemental  filters,  the  differences  are  negligible  and  the  results  of  Sections  3.3.1  and  3.7.1 
apply  to  this  application.  While  Section  3.7.1  shows  that  the  transformation  matrix  Tj  (ti) 
can  be  eliminated  from  the  hypothesis  conditional  probability  calculation,  it  will  be  used 
in  the  results  presented  in  this  section.  Thus,  the  effect  of  neglecting  the  differences  in 
Tj  (ti)  will  be  seen.  The  analysis  could  have  been  repeated  implementing  the  assumption 
that  Tj  (ti)  does  not  vary  among  the  elemental  filters  and  removing  it  from  hypothesis  con¬ 
ditional  probability  calculation.  However,  the  same  result  (in  the  next  paragraph)  would 
have  been  achieved.  Therefore,  Tj  (ti)  for  all  elemental  filters  was  used  in  the  analysis  that 
follows. 

To  demonstrate  that  any  7  value  results  in  equivalent  hypothesis  conditional  proba¬ 
bilities,  a  test  was  performed.  GPS  double  difference  widelane  carrier-phase  measurements 
were  used  in  GRMMAEs,  with  7  =  0  and  7  =  1  (i.e.,  using  post-fit  and  traditional  residu¬ 
als,  respectively).  The  double  difference  measurements  were  taken  over  a  69  km  baseline. 
In  both  cases,  100  elemental  filters  were  built  upon  100  hypothesized  ambiguity  vectors 
generated  using  the  Fast  Ambiguity  Search  Filter  (FASF)  routine  [7].  The  hypothesis  con¬ 
ditional  probability  for  each  elemental  filter  was  calculated  in  each  of  the  GRMMAFs.  A 
difference  in  probabilities  was  then  calculated  by  subtracting  the  hypothesis  conditional 
probabilities  calculated  in  the  traditional  MMAE  (GRMMAE  with  7  =  1)  from  those  in 
the  post-fit  modified  MMAE  (GRMMAE  with  7  =  0).  The  probability  difference  between 
the  two  GRMMAEs  was  on  the  order  of  10“^^  and  was  likely  due  to  machine  precision. 
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rounding  errors,  and  the  slight  differences  in  measurement  matrices  previously  mentioned. 
Since  the  probability  for  the  correct  hypothesis  is  on  the  order  of  1,  a  difference  of  10“^^  is 
essentially  zero.  Thus,  this  demonstrates  that  the  GRMMAE  with  7  =  1  is  equivalent  to 
the  GRMMAE  with  7  =  0  for  this  application,  as  predicted  by  the  theory  in  Section  3.7.1. 


4-3  Beta  Dominance 

This  dissertation  repeatedly  states  that  the  GRMMAE  is  susceptible  to  /3-dominance, 
described  in  Section  2.4.1,  when  the  transformation  matrix  T j  {ti)  varies  among  the  ele¬ 
mental  filters.  Derivation  of  this  fact  in  Section  3.6  is  one  of  the  contributions  of  this 
research.  This  section  will  demonstrate  that  fact  by  showing  an  example  in  which  the 
GRMMAE  with  7  =  1  (i.e.,  a  traditional  MMAE)  functions  properly  while  a  GRMMAE 
with  7  =  0  (i.e.,  a  post-fit  modified  MMAE)  experiences  /3-dominance. 

Recall  from  Section  3.6  that  the  multivariate  Gaussian  density  function  used  to  cal¬ 
culate  the  hypothesis  conditional  probabilities  in  the  GRMMAE  is  given  by 


1 

(27r)^abs(|Tj(3i)  \)\Aj  (t”)  |5 


exp|-^rj  (tj  A/  {t-  )| 


(4.16) 


This  density  is  composed  of  a  “beta  term”  and  an  “exponential  term,”  respectively,  given 
by 


1 

(27r)2  abs(|Tj(ti)|)|Aj  (t“)  |2 

)  ^7^  (7“ )  )  }  (4-17) 

In  Equation  (4.17),  assuming  linear  system  models  are  used,  only  the  exponential  term 
depends  upon  the  measurements.  The  beta  term  is  based  solely  on  the  system  models.  As 
such,  it  is  completely  pre-computable  before  any  measurements  are  taken.  Therefore,  to 
ensure  the  GRMMAE  converges  to  the  correct  elemental  filter,  the  hypothesis  conditional 
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probability  calculation  must  be  driven  by  the  exponential  term  and  not  the  beta  term.  In 
some  applications,  the  system  models  are  such  that  the  beta  terms  have  a  wide  separation 
between  them,  and  the  residuals  are  similar  enough  to  cause  the  exponential  terms  to 
have  similar  magnitudes.  When  this  happens,  the  beta  terms  dominate  the  hypothesis 
conditional  probability  calculation,  and  incorrect  results  may  be  achieved. 

An  example  of  /3-dominance  is  created  here  to  illustrate  the  GRMMAE’s  susceptibil¬ 
ity  to  this  phenomenon  and  to  explain  /3-dominance  further.  The  truth  model  of  Equations 
(4.3)  through  (4.12)  was  used  with  a  slight  change  to  Equation  (4.10).  Eor  this  example, 
the  true  noise  covariance  matrices  were 


0.0454 

0.0838 

0.0646 

R  (ti)  = 

0.0025 

0 

0.0838 

0.1548 

0.1193 

0 

0.0025 

0.0646 

0.1193 

0.0920 

(4.18) 


The  measurement  noise  covariance  in  Equation  (4.18)  was  the  same  as  that  in  High  Q^^/Low 
R  scenario  of  Equation  (4.10),  while  the  equivalent  discrete-time  noise  covariance  was  2.5 
times  that  in  the  High  Q^^/Low  R  scenario  of  Equation  (4.10).  The  change  to  the  equivalent 
discrete-time  noise  covariance  was  made  to  facilitate  creation  of  /3-dominance. 

The  application  used  for  this  example  is  the  unknown  application,  similar  to 
Section  4.1.3.  Three  elemental  filters  were  built  based  upon  three  hypotheses  for  the 
equivalent  discrete-time  noise 


2Qd  (ti)  2.5Qd{ti)  SQditi) 

where  is  the  matrix  given  in  the  High  Q^j/Low  R  scenario  of  Equation  (4.10).  Since 
there  is  only  a  small  change  between  these  three  hypotheses,  it  is  expected  that  the  GR- 
MMAE  will  converge  to  the  correct  solution  relatively  slowly.  R  is  also  expected  that,  if 
there  is  a  large  separation  between  the  beta  terms  in  the  elemental  filters,  then  the  GR- 
MMAE  will  experience  /3-dominance.  Here,  “a  large  separation”  is  most  easily  explained 
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graphically.  A  GRMMAE  was  simulated  100  times  with  7  =  1  (i.e.,  a  traditional  MMAE), 
and  each  simulation  consisted  of  500  measurement  epochs.  Eigure  4.12  is  a  plot  of  the 
sample  averages  of  the  exponential  terms  and  the  beta  terms  for  these  simulations,  and 
Eigure  4.13  gives  an  expanded  plot  of  the  exponential  terms.  Since  the  beta  terms  do  not 
depend  upon  the  measurements,  they  are  the  same  in  all  100  simulations. 

In  Eigure  4.12  the  exponential  terms  all  have  similar  magnitudes.  Recall  from  Section 
2.4.3,  that  it  is  desirable  for  to  be  on  the  order  of  the  measure¬ 

ment  dimension  (2  in  this  example)  [46].  Thus,  the  exponential  term  for  the  correct  filter 
should  be  on  the  order  of  exp(— 1)  =  0.368,  which  it  is  in  this  example.  However,  the 
exponential  terms  from  the  other  elemental  filters  are  on  that  same  order.  As  a  result,  the 
GRMMAE  with  7  =  1  converges  rather  slowly  to  the  correct  hypothesis  (as  will  be  shown 
in  Eigure  4.17).  Also  note  that  the  beta  terms  in  Eigure  4.12  are  relatively  small  and  closely 
spaced.  Thus,  it  is  expected  that  this  GRMMAE  (with  7  =  1)  will  not  be  susceptible  to 
/3-dominance.  To  continue  with  this  illustration,  Eigure  4.14  shows  a  plot  of  the  product  of 
the  beta  terms  and  the  sample  average  exponential  terms  at  each  measurement  epoch.  In 
Eigure  4.14,  part  (a)  shows  the  full  plot  while  part  (b)  expands  on  the  early  portion  of  the 
run  to  show  more  detail.  Erom  Eigure  4.14,  the  products  of  the  beta  and  sample  average 
exponential  terms  have  similar  magnitudes  for  all  three  elemental  filters.  Therefore,  the 
GRMMAE  hypothesis  conditional  probability  calculation  is  not  dominated  by  the  beta 
terms.  Rather,  the  GRMMAE  is  able  to  use  both  residuals  (i.e.,  the  exponential  terms) 
and  the  prior  probabilities  in  the  hypothesis  conditional  probability  calculation. 

In  contrast,  the  same  100  measurement  sets  were  used  to  repeat  this  example  for  a 
GRMMAE  with  7  =  0  (i.e.,  a  post-fit  modified  MMAE).  In  addition  to  using  the  same 
measurements,  the  same  GRMMAE  was  used  with  the  exception  that  the  7  design  pa¬ 
rameter  was  changed  from  1  to  0.  Thus,  a  direct  comparison  between  GRMMAEs  with 
7  =  1  and  7  =  0  was  made.  Eigure  4.15  shows  a  plot  of  the  beta  and  the  sample  average 
exponential  terms  for  the  GRMMAE  with  7  =  0.  Gomparing  Eigures  4.12  and  4.15  shows 
that  the  sample  average  exponential  terms  were  unchanged  by  the  7  design  parameter, 
as  derived  in  Section  3.6.  This  confirms  that  the  only  dependence  of  the  GRMMAE  hy¬ 
pothesis  conditional  probability  calculation  on  the  7  design  parameter  is  caused  by  the 
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Figure  4.12  Beta  and  Sample  Average  Exponential  Terms  for  GRMMAE  with  7 
(Traditional  MMAE) 
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Exponential  Term  versus  Measurement  Number  (expanded) 
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Figure  4.13  Sample  Average  Exponential  Terms  (expanded)  for  GRMMAE  with  7  =  1 

(Traditional  MMAE) 

beta  terms,  demonstrating  one  of  the  key  contributions  of  this  dissertation.  Notice  that, 
in  this  example,  the  beta  terms  for  the  GRMMAE  with  7  =  0  are  two  orders  of  magni¬ 
tude  larger  than  the  beta  terms  for  the  GRMMAE  with  7  =  1.  This  was  caused  by  the 
determinant  of  the  transformation  matrix  in  the  denominator  of  the  beta  term.  When 
7=1,  |Tj  (tj)  I  =  1  Vj  and  it  has  no  effect  on  the  beta  term.  However,  when  7  7^  1,  the 
determinant  has  an  effect  on  the  beta  terms.  In  this  example,  7  =  0  causes  abs(|Tj  (tj)  |) 
to  increase  greatly  the  magnitude  of  the  beta  terms^. 

In  addition  to  increasing  the  magnitude  of  the  beta  terms,  the  separation  between 
the  beta  terms  was  also  increased  when  7  =  0.  In  this  example,  the  increase  in  separation 
was  enough  to  cause  /3-dominance.  Figure  4.16,  a  plot  of  the  product  of  the  beta  and 
sample  average  exponential  terms,  demonstrates  this  /3-dominance.  Unlike  in  Figure  4.14, 
the  product  of  the  beta  and  the  sample  average  exponential  terms  shown  in  Figure  4.16 
has  sufficient  separation  between  the  hypotheses  to  cause  the  beta  terms  to  dominate  over 

^It  is  not  correct  to  say  that  abs(|T  (7)  |)  amplifies  the  beta  terms  because  the  shape  of  the  beta  plots 
is  changed  in  addition  to  the  magnitude. 
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(a)  Beta  Terms  (b)  Exponential  Terms 

Figure  4.15  Beta  and  Sample  Average  Exponential  Terms  for  GRMMAE  with  7  =  0 
(Traditional  MMAE) 

the  exponential  terms.  At  every  measurement  epoch  in  Figure  4.16  (b),  the  product  of  the 
beta  and  sample  average  exponential  terms  for  elemental  filter  3  was  larger  than  the  same 
product  for  the  other  filters.  This  was  caused  by  the  dominant  influence  of  the  beta  term 
for  elemental  filter  3.  The  root  cause  of  this  was  the  separation  between  the  beta  terms 
created  by  abs(|Tj  (U)  |). 

As  final  evidence  that  the  GRMMAE  with  7  =  0  suffered  from  /3-dominance  while 
the  one  with  7  =  1  did  not,  the  sample  mean  hypothesis  conditional  probabilities  from 
the  100  Monte  Garlo  runs  are  shown  in  Figure  4.17.  Part  (a)  of  the  figure  shows  the 
sample  average  hypothesis  conditional  probabilities  for  the  GRMMAE  with  7  =  1  (i.e.,  a 
traditional  MMAE).  In  this  figure,  the  probability  for  the  correct  elemental  filter  slowly 
converged  toward  1.  It  is  reasonable  to  assume  that  this  probability  would  have  converged 
to  1  if  the  simulation  had  been  allowed  to  continue  longer.  Part  (b)  of  the  figure,  on 
the  other  hand,  shows  the  hypothesis  conditional  probabilities  for  the  GRMMAE  with 
7  =  0  (i.e.,  a  post-fit  modified  MMAE).  In  that  figure,  the  probability  for  elemental 
filter  3,  an  incorrect  hypothesis,  very  quickly  converged  to  1,  while  the  probability  for  the 
correct  hypothesis  converged  to  0.  This  occurred  because  the  products  of  the  beta  and 
sample  average  exponential  terms,  shown  in  Figure  4.16,  for  the  hypotheses  were  separated 
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Figure  4.16 
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Beta  Term  Times  Exponential  Term 


(b)  Beta  Term  Times  Exponential  Term  (expanded) 


Product  of  Beta  and  Sample  Average  Exponential  Terms  for  GRMMAE 
with  7  =  0  (Post-Eit  Modified  MMAE) 
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Probability  Probability 


Figure  4.17 


Hypothesis  Conditional  Probabilities 


(a)  Probabilities  for  GRMMAE  with  7=1 


Hypothesis  Conditional  Probabilities 


(b)  Probabilities  for  GRMMAE  with  7  =  0 


Hypothesis  Conditional  Probability  Comparison  for  GRMMAEs  with  7 
and  7  =  0 
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sufficiently  to  allow  the  beta  terms  to  dominate  the  hypothesis  conditional  probability 
calculation.  In  this  example,  with  7  =  0,  elemental  filter  3  had  the  largest  beta  term,  as 
shown  in  Figure  4.12.  Therefore,  the  /3-dominance  caused  the  GRMMAE  with  7  =  0  to 
converge  to  elemental  filter  3,  despite  the  fact  that  this  was  an  incorrect  hypothesis. 

4-4  Effect  of  'j  on  Beta  Terms 

The  previous  section  demonstrated  that  the  design  parameter  7  changes  the  magni¬ 
tude  and  the  shape  of  the  (3*  (tj)  plots.  A  mathematical  derivation  of  the  effects  of  7  on 
P*  (ti)  does  not  support  the  primary  goal  of  this  research  and,  therefore,  is  not  pursued. 
However,  certain  properties  were  observed  and  are  demonstrated  here.  These  observations 
are  given  without  derivation  or  proof  in  support  of  a  recommendation  for  future  research 
given  in  Chapter  V. 

First,  it  was  noted  in  the  previous  section  that  7  changes  the  shape  of  the  /3*  {tp 
plots.  This  occurs  only  in  the  transient  portion  of  the  curve  before  the  elemental  filters 
reach  steady  state.  For  most  7  values,  /3t  (U)  increases  to  a  steady  state  value  during  the 
transient  portion  of  the  curve.  This  behavior  is  seen  in  Figure  4.12  (a).  This  shape  makes 
sense  because  the  residual  covariance  Aj  (t~),  which  appears  in  the  denominator  of  /3*  {tp, 
typically  decreases  as  the  elemental  filter  reaches  steady  state.  However,  at  some  point 
as  7  approaches  zero,  the  curve  flips  such  that  /3*  {f)  decreases  to  a  steady  state  value 
during  the  transient  portion  of  the  curve,  as  seen  in  Figure  4.15  (a).  This  flip  appears  to 
occur  when  abs(|Tj  {tp  |)  begins  to  dominate  over  |Aj  in  the  denominator  of  the 

beta  terms.  Analogous  to  /3-dominance,  this  effect  could  be  termed  “T-dominance.”  This 
T-dominance  effect  appears  to  occur  when  7  is  decreased  below  a  threshold  value  which 
causes  abs(|Tj  (ti)  |)  to  become  small  enough  to  dominate  over  |A  .  A)  t'"  in  the  beta 
term.  Confirmation  of  this  observation  and  theoretical  determination  of  the  threshold  7  are 
not  pursued  here,  since  they  do  not  support  the  primary  goal  of  this  research.  However,  if  a 
theoretical  threshold  for  7  is  determined,  it  may  be  possible  to  use  the  CRMMAE  to  correct 
for  /3-dominance  rather  than  using  the  traditional  method  of  /3  stripping.  Alternately,  this 
flip  in  the  /3*  (f)  curve  may  be  related  to  the  unexpected  shape  of  Eigure  4.2,  in  which 
the  generalized  residual  covariance  increased  to  a  steady  state  value  when  7  =  0  (i.e.. 
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post-fit  residual)  was  used.  Again,  there  may  be  a  threshold  7  which  caused  that  behavior. 
Finding  this  7  threshold  would  support  development  of  a  method  to  use  7  to  correct  for 
/3-dominance. 

Secondly,  it  was  observed  that,  as  7  becomes  large  in  absolute  value,  the  beta  terms 
generally  become  small.  To  demonstrate  one  example  of  this,  the  analysis  from  the  previous 
section  was  re-run  for  7  equal  100,  10,  -10,  and  -100,  and  the  beta  terms  were  plotted  in 
Figure  4.18.  Since  the  beta  terms  do  not  vary  from  one  Monte  Carlo  run  to  the  next, 
there  was  no  need  to  run  multiple  simulations  for  this  example.  The  first  observation  from 
Figure  4.18  is  that  the  beta  plots  were  so  similar  that  they  appear  the  same  for  numbers  of 
equal  magnitude  but  opposite  signs  (i.e.,  100  and  -100  produce  nearly  identical  plots).  The 
second  and  more  important  observation  is  that  larger  (magnitude)  7  values  led  to  smaller 
(3*  {ti)  values.  This  inverse  relationship  is  also  demonstrated  in  Figure  4.18. 

From  Figure  4.18,  it  appears  that  increasing  7  by  an  order  of  magnitude  causes  the 
beta  terms  to  decrease  by  an  order  of  magnitude.  This  observation  is  further  supported 
by  comparing  Figures  4.18  and  4.12.  However,  Figure  4.15  clearly  shows  that  relationship 
does  not  hold  for  all  7.  From  Figures  4.12  and  4.15,  decreasing  7  from  1  to  0  caused  the 
beta  terms  to  increase  by  2  orders  of  magnitude.  Thus,  some  type  of  singular  behavior 
is  observed  as  7  ^  0.  Again,  this  observation  does  not  support  the  primary  goal  of  this 
research.  Therefore,  a  mathematical  description  of  this  behavior  is  not  pursued  in  this 
dissertation.  However,  a  mathematical  description  of  the  relationship  between  7  and  the 
beta  terms  is  probably  the  key  to  using  the  GRMMAE,  with  an  appropriate  7,  to  correct 
/3-dominance,  rather  than  using  the  ad  hoc  /3  stripping. 

4-5  Scalar  Penalty  Increase 

Section  2.4.5  mentioned  scalar  penalty  increase  as  an  alternative  to  the  GRMMAE 
to  cause  faster  MMAE  convergence  to  the  best  parameter  [31].  Although  scalar  penalty 
increase  is  not  a  research  contribution  of  this  dissertation,  an  example  is  presented  here, 
since  it  is  an  alternative  to  the  GRMMAE.  Monte  Garlo  simulations  of  a  system  with  an 
unknown  measurement  bias  were  used  to  demonstrate  scalar  penalty  increase.  With  the 
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Figure  4.18  Beta  Terms  for  GRMMAEs  with  7  =  ±100  and  7  =  ±10 
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exception  of  the  noise  covariance  matrices,  the  truth  model  from  Section  4.1.1  was  used 


for  these  simulations.  The  noise  covariance  matrices  for  these  simulations  were 


0.0454 

0.0838 

0.0646 

_ 

_ 

R  (ti)  = 

1.5625 

0 

0.0838 

0.1548 

0.1193 

X  10“^ 

0 

1.5625 

0.0646 

0.1193 

0.0920 

(4.19) 


As  before,  three  elemental  filters  were  used.  The  three  hypotheses  for  the  measure¬ 
ment  bias  were 


0 

-1 

1 

n2  = 

ns  = 

0 

0 

0 

(4.20) 


The  true  bias  was  ni.  These  hypotheses  were  changed  from  those  in  Section  4.1.2  to  create 
a  more  challenging  scenario  for  the  MMAE.  By  making  them  very  similar,  the  residuals 
from  the  three  elemental  filters  will  also  be  similar.  This  causes  the  MMAE  to  converge 
to  the  correct  hypothesis  more  slowly,  allowing  a  better  opportunity  to  demonstrate  the 
effect  of  scalar  penalty  increase. 


The  system  was  simulated  1000  times  each  for  a  traditional  MMAE  (i.e.,  GRMMAE 
with  7=1)  with  and  without  scalar  penalty  increase.  Figure  4.19  (a)  shows  the  sample 
average  hypothesis  conditional  probabilities  calculated  for  the  traditional  MMAE  without 
scalar  penalty  increase.  Figure  4.19  (b)  shows  the  hypothesis  conditional  probability  for 
the  traditional  MMAE  with  scalar  penalty  increase.  In  the  MMAE  with  scalar  penalty 
increase,  the  exponential  term  in  the  Gaussian  probability  densities  of  Equation  (2.41)  was 
changed  from  exp{  (t” )  (t” )  (t”)}  to  exp{  -2r^  (t” )  (t” )  (t”)}. 

Figure  4.19  demonstrates  that  the  MMAE  converged  to  the  correct  measurement 
bias  vector  regardless  of  whether  scalar  penalty  increase  was  used  or  not.  However,  the 
MMAE  with  scalar  penalty  increase  converged  more  quickly  than  the  MMAE  without 
scalar  penalty  increase.  This  simple  example  demonstrates  that  scalar  penalty  increase 
is  useful  to  cause  faster  convergence,  the  goal  the  GRMMAE  was  originally  proposed  to 
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Figure  4.19 


Hypothesis  Conditional  Probability 


(a)  Traditional  MMAE  Without  Scalar  Penalty  Increase  (a  =  1/2) 


Hypothesis  Conditional  Probability 


(b)  Traditional  MMAE  With  Scalar  Penalty  Increase  {a  =  2) 


Hypothesis  Conditional  Probability  for  a  Traditional  MMAE  With  and 
Without  Scalar  Penalty  Increase 
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accomplish.  The  drawback  to  scalar  penalty  increase,  however,  is  that  it  can  lead  to  false 
alarms,  as  discussed  in  Section  2.4.5.  Additionally,  the  “increased  scalar”  value  must  be 
chosen  in  an  ad  hoc  manner.  There  is  no  method  for  choosing  an  optimal  scalar  value. 

4-6  Chapter  Summary 

This  chapter  provided  simulations  to  verify  the  theory  derived  in  Chapter  III.  First, 
the  generalized  residual  vector  and  covariance  matrix  equations  from  Section  3.2  were  ver- 
ihed  using  a  series  of  simple  simulations.  Next,  since  GPS  carrier-phase  integer  ambiguity 
resolution  is  the  application  for  which  the  GRMMAE  was  originally  designed,  double  dif¬ 
ference  carrier-phase  measurements  were  used  to  demonstrate  that  the  GRMMAE,  when 
applied  to  a  measurement  bias  problem  (same  T  (tj)  in  all  elemental  hlters),  is  independent 
of  the  7  design  parameter.  Then,  the  unknown  application  (differing  Tj  (U)  in  elemen¬ 
tal  hlters)  was  used  to  demonstrate  that  the  MMAE  is  susceptible  to  the  /3-dominance 
effect.  This  example  demonstrated  how  the  7  parameter  created  the  /3-dominance  effect. 
Next,  observations  concerning  the  effect  of  7  on  the  /3t  (tj)  terms  were  described.  These 
observations  are  likely  to  serve  as  the  base  point  for  recommended  future  research  into  use 
of  7  to  correct  for  /3-dominance.  Finally,  scalar  penalty  increase  was  demonstrated  as  an 
alternative  to  the  GRMMAE.  With  the  exception  of  Sections  4.4  and  4.5,  all  sections  in 
this  chapter  demonstrate  scientihc  contributions  of  this  dissertation. 
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V.  Conclusions  and  Recommendations 


This  dissertation  has  developed  the  mathematical  foundation  for  using  generalized  resid¬ 
uals  to  calculate  hypothesis  conditional  probabilities  in  a  Multiple  Model  Adaptive  Es¬ 
timator  (MMAE).  The  developed  estimator  is  termed  a  Generalized  Residual  Multiple 
Model  Adaptive  Estimator  (GRMMAE).  This  GRMMAE  was  inspired  by  the  use  of  post- 
fit  residuals  in  a  modified  MMAE  for  GPS  carrier-phase  integer  ambiguity  resolution. 
Previous  research  demonstrated  that  the  modified  MMAE  can  successfully  resolve  GPS 
carrier-phase  integer  ambiguities  [17,  18].  However,  that  research  did  not  provide  a  math¬ 
ematical  derivation  of  the  modified  MMAE  or  the  effects  of  the  post-fit  residuals  on  the 
MMAE’s  performance.  Thus,  in  an  effort  to  expand  post-fit  residual  use  in  MMAE  ap¬ 
plications,  this  research  provided  a  derivation  of  the  GRMMAE.  Through  this  derivation, 
the  effect  of  post-fit  residuals  on  the  GRMMAE’s  performance  has  been  determined.  This 
final  chapter  provides  a  chapter- by-chapter  summary  of  the  dissertation.  Section  5.1  lists 
the  research  contributions  made  in  this  dissertation.  Then,  the  chapter  concludes  with  a 
recommendation  for  future  research. 

Ghapter  I  of  this  dissertation  presented  an  introduction  to  the  problem  of  estimating 
system  states  in  the  presence  of  an  unknown  system  parameter.  It  discussed  common  prob¬ 
lems  in  which  system  parameters  must  be  estimated  in  addition  to  estimating  the  system 
state,  and  gave  the  motivation  for  using  the  MMAE  for  parameter  and  state  estimation. 
Einally,  it  discussed  the  modified  MMAE,  which  uses  post-fit  residuals,  as  the  motivation 
for  this  research. 

Ghapter  II  presented  the  background  theory  necessary  to  derive  the  GRMMAE. 
It  presented  derivations  of  both  the  linear  and  nonlinear  (extended)  Kalman  filters.  It 
gave  a  derivation  of  the  traditional  MMAE  and  described  several  common  performance 
enhancements  to  the  MMAE.  Since  this  dissertation  was  motivated  by  a  GPS  application, 
it  also  presented  a  brief  introduction  to  the  GPS.  It  discussed  precise  GPS  and  the  need 
to  resolve  carrier-phase  integer  ambiguities  in  order  to  estimate  precise  relative  positions 
with  GPS.  Ghapter  II  concluded  by  describing  the  modified  MMAE  for  GPS  carrier-phase 
integer  ambiguity  resolution.  The  modified  MMAE  is  the  technique  which  motivated  this 
dissertation. 
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Chapter  III  presented  the  research  contributions  of  this  dissertation.  This  chapter 
provided  the  mathematical  derivation  of  the  GRMMAE.  This  GRMMAE  is  a  new  Multiple 
Model  Adaptive  Estimator  designed  such  that,  by  varying  a  single,  scalar  design  parame¬ 
ter,  the  GRMMAE  can  be  equivalent  to  a  traditional  MMAE,  a  modified  MMAE,  or  any 
other  MMAE  based  on  a  linear  combination  of  traditional  and  post-fit  Kalman  filter  resid¬ 
uals.  Chapter  III  began  with  the  definition  of  the  new  generalized  residual.  Next,  it  gave  a 
derivation  of  the  properties  of  the  generalized  residual  vector  from  the  elemental  Kalman 
filter  based  on  the  correct  system  parameter.  After  that,  the  properties  for  the  generalized 
residual  vectors  from  incorrect  elemental  filters  were  derived  for  several  common  MMAE 
applications.  Chapter  III  showed  that  the  generalized  Kalman  filter  residual  is  related  to 
the  traditional  Kalman  filter  residual  through  a  simple  transformation  matrix  T  (tj).  The 
effect  of  this  transformation  matrix  was  described  in  Chapter  III,  and  it  gave  a  derivation  of 
the  GRMMAE  hypothesis  conditional  probability  formula  based  on  the  generalized  resid¬ 
ual.  The  relationship  between  the  traditional  MMAE  hypothesis  conditional  probability 
formula  and  the  GRMMAE  hypothesis  conditional  probability  formula  was  presented.  This 
relationship  demonstrated  that  the  GRMMAE  design  parameter  7  effects  only  the  (3*  {ti) 
term  in  the  GRMMAE  hypothesis  conditional  probability  calculation.  This  finding  is  the 
most  significant  finding  in  this  dissertation,  since  it  means  that  the  GRMMAE  is  artificially 
susceptible  to  /3-dominance.  Chapter  III  concluded  by  further  developing  the  generalized 
hypothesis  conditional  probability  formula  under  two  conditions:  1)  the  transformation 
matrix  T  {ti)  is  the  same  in  all  elemental  filters;  2)  the  transformation  matrix  T  {ti)  is  not 
the  same  in  all  elemental  filters.  It  showed  that,  under  condition  1),  the  design  parameter 
7  is  cancelled  from  the  hypothesis  conditional  probability  calculation.  Thus,  under  con¬ 
dition  1),  the  GRMMAE  produces  results  equivalent  to  a  traditional  MMAE  regardless 
of  the  value  of  7.  In  that  case,  the  GRMMAE  is  not  susceptible  to  /3-dominance.  Under 
condition  2),  on  the  other  hand,  the  effect  of  7  on  the  hypothesis  conditional  probability 
calculation  does  not  cancel.  In  that  case,  the  GRMMAE  is  susceptible  to  /3-dominance. 

Chapter  IV  presented  simulation  results  to  demonstrate  the  theory  derived  in  Chap¬ 
ter  III.  It  began  with  a  series  of  simulations  to  verify  the  equations  for  calculating  the 
generalized  residual  vector  and  its  covariance  matrix  from  the  traditional  residual  vector. 


5-2 


Next,  the  GRMMAE  was  used  to  demonstrate  that  applications  with  the  same  T  (tj)  in 
all  elemental  filters  are  independent  of  7.  Since  the  GRMMAE  derivation  was  motivated 
by  the  modified  MMAE  for  GPS  carrier-phase  integer  ambiguity  resolution,  GPS  carrier- 
phase  measurements  were  used  for  this  example.  Next,  the  susceptibility  of  the  GRMMAE 
to  /3-dominance  was  demonstrated.  Then,  the  effect  of  7  on  the  /3t  (t*)  terms  was  discussed. 
The  observations  about  this  effect  represent  a  starting  point  for  the  recommended  future 
research  presented  later  in  this  chapter. 

5.1  Research  Contributions 

The  development  of  the  Generalized  Residual  Multiple  Model  Adaptive  Estimator  is 
the  primary  contribution  of  this  dissertation.  The  heart  of  this  development  is  the  deriva¬ 
tion  of  the  hypothesis  conditional  probability  calculation  based  on  generalized  residuals. 
This  development  provides  a  simple,  general  design  in  which  a  single,  scalar  design  parame¬ 
ter  7  is  used  to  specify  a  traditional  MMAE,  a  modified  MMAE,  or  any  other  MMAE  using 
a  linear  combination  of  traditional  and  post-fit  residuals.  Several  supporting  contributions 
are  developed  in  this  dissertation  as  part  of  the  GRMMAE  foundation.  These  support¬ 
ing  contributions  are  derived  in  Ghapter  III  and  verified  in  Ghapter  IV.  The  supporting 
contributions  are: 

•  Derived  generalized  residual  properties 

•  Derived  and  described  transformation  T  (tj)  from  traditional  to  generalized  residual 

•  Derived  hypothesis  conditional  probability  formula  using  generalized  residual 

•  Derived  relationship  between  traditional  and  generalized  hypothesis  conditional  prob¬ 
ability  formulas 

•  Derived  and  demonstrated  equivalence  between  MMAE  and  GRMMAE  when  ele¬ 
mental  filters  have  same  T  (tj) 

•  Derived  formula  demonstrating  GRMMAE  susceptibility  to  /3-dominance 

•  Demonstrated  how  7  7^  1  causes  /3-dominance 

Each  of  these  contributions  are  summarized  in  the  sections  that  follow. 
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5.1.1  Derived  Generalized  Residual  Properties.  Generalized  residual  properties 
were  derived  in  Sections  3.2  and  3.3.  Section  3.2  demonstrated  that  the  generalized  residu¬ 
als  from  the  correct  elemental  filter  are  a  linear  function  of  traditional  residuals  from  that 
same  elemental  filter.  Since  the  traditional  residuals  are  white,  Gaussian,  and  zero-mean 
[29],  the  generalized  residuals  from  the  correct  filter  are  also  white,  Gaussian,  and  zero- 
mean.  Additionally,  given  that  is  the  covariance  of  the  traditional  residuals  from 

the  correct  elemental  filter  [29,  30],  the  covariance  of  the  generalized  residuals  from  the 
correct  elemental  filter  is  given  in  Equation  (3.4),  repeated  here  for  convenience 


AliU)  =  T,{U)Ak{t-)Tl{U)  (3.4) 

where  (fj)  =  [l  —  (1  —  7)  (t*)  (tj)]  is  the  linear  relationship  between  the  tra¬ 

ditional  and  generalized  residuals  and  7  is  defined  in  Equation  (3.1). 

In  Section  3.3,  the  residual  vector  mean  and  covariance  for  an  elemental  filter  with 
an  unknown  measurement  bias  were  shown  to  be 

m^.  (ti)  =  T  (ti)  [H  (ti)  (xfc  {t~ )  -  Xj  {t~))  +  Arifcj]  (3.11) 

A*{U)  =T{U)Ak{t-)T^{ti)  (3.13) 

In  this  case,  it  was  demonstrated  that  H  (ti)  and  K  (ti)  were  the  same  in  all  elemental  filters. 
Thus,  T  (ti)  was  the  same  in  all  elemental  filters  and  the  subscript  identifying  the  elemental 
filter  was  removed.  Equations  (3.11)  and  (3.13)  were  shown  to  apply  to  any  application  in 
which  the  transformation  matrix  T  (ti)  does  not  vary  among  the  elemental  filters.  Eor  an 
unknown  (ti),  F  (t)  (and  thus  $  (tj+i,  tj)),  R  (U),  or  any  other  application  for  which  the 
transformation  matrix  T  (ti)  varies  among  the  elemental  filters,  the  generalized  residual 
mean  and  covariance  for  an  incorrect  filter  were  shown  to  be 

mr-  (ti)  =  Tj  (ti)  H  (ti)  [Sfc  (t- )  -  %  (t- )]  (3.15) 

A*(U)  =T,(U)Ak(t-)Tj(U) 
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These  residual  properties  support  the  main  contribution  of  providing  a  mathematical  foun¬ 
dation  for  the  use  of  generalized  residuals  in  an  MMAE.  Similar  derivations  of  traditional 
residual  properties  for  a  traditional  MMAE  are  found  in  [16] . 

5.1.2  Derived  and  Deseribed  Transformation  From  Traditional  to  Generalized  Resid¬ 
ual  T(tj).  The  transformation  from  traditional  to  generalized  residuals  is  given  in 
Equation  (3.2).  It  is 

v*{U)  =  T{U)v{t-)  (3.2) 

Section  3.4  demonstrates  that  the  effect  of  this  transformation  is  to  scale  and  rotate  the 
residual  vector  in  the  measurement  space.  Additionally,  it  was  stated  that,  in  general, 
the  matrix  T  (tj)  is  not  idempotent,  so  it  is  a  transformation  and  not  a  projection.  These 
contributions  support  the  primary  contribution  by  deriving  and  characterizing  the  effect 
of  generalized  residuals  on  the  GRMMAE  hypothesis  conditional  probability  calculation. 

5.1.3  Derived  Hypothesis  Conditional  Probability  Formula  Using  Generalized  Resid¬ 
ual.  The  derivation  of  the  hypothesis  conditional  probability  formula  using  generalized 
residuals  is  given  in  Section  3.5.  This  derivation  is  the  heart  of  this  dissertation.  Previous 
research  implemented  the  MMAE  with  post-fit  residuals  [17,  18].  However,  that  research 
did  not  mathematically  derive  the  hypothesis  conditional  probability  formula,  nor  did  it 
mathematically  characterize  the  effect  of  using  non-traditional  residuals  in  the  MMAE. 
In  Section  3.5,  this  dissertation  provided  a  mathematical  derivation  of  the  GRMMAE  hy¬ 
pothesis  conditional  probability  calculation.  It  was  shown  that,  by  setting  7  =  0,  the 
GRMMAE  is  equivalent  to  the  modified  MMAE  using  post-fit  residuals.  However,  the 
GRMMAE  is  more  general  than  the  modified  MMAE.  By  setting  7  =  1,  the  GRMMAE  is 
equivalent  to  a  traditional  MMAE.  Additionally,  by  setting  7  equal  to  any  other  real,  scalar 
number,  the  GRMMAE  implements  an  infinite  number  of  new  MMAE’s  not  encountered 
in  scientific  literature  to  date.  Thus,  the  GRMMAE  derivation  provides  a  mathematical 
foundation  for  the  modified  (i.e.,  post-fit  residual)  MMAE  which  also  applies  to  a  more 
general  class  of  MMAE’s.  This  mathematical  foundation  is  not  found  in  scientific  literature 
to  date  and  is  the  primary  contribution  of  this  dissertation. 
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5.1.4  Derived  Relationship  Between  Traditional  and  Generalized  Hypothesis  Condi¬ 

tional  Probability  Formulas.  Using  the  GRMMAE  mathematical  foundation,  the  effect  of 
generalized  residuals  in  the  MMAE  is  derived  in  Section  3.6.  This  derivation  demonstrates 
that  using  7  7^  1  affects  only  the  (3*  {ti)  terms  in  the  GRMMAE  hypothesis  conditional 
probability  calculation.  The  specific  choice  of  7  will  increase  or  decrease  the  magnitudes  of 
the  l3*  {ti)  terms.  Thus,  7  is  similar  to  a  tuning  parameter.  It  is  used  to  “tune”  the  relative 
weight  of  the  (3*  (ti)  terms  compared  to  the  exp  |  — l/2rj  rj  (t7)|  terms  in 

the  hypothesis  conditional  probability  calculation.  However,  unlike  the  exponential  terms, 
the  l3j  {ti)  terms  are  based  solely  on  the  elemental  filter  models  and  are  unaffected  by  the 
real-world  measurements  (assuming  the  elemental  filter  models  are  linear).  Therefore,  us¬ 
ing  7  to  artificially  increase  the  f3j  {U)  terms  is  ill-advised,  as  it  may  lead  to  /3-dominance 
and  incorrect  results  in  the  GRMMAE.  This  finding  is  a  key  contribution,  as  it  points  out 
a  potentially  fatal  flaw  if  other  researchers  attempt  to  use  post-fit  residuals  (or  other  gen¬ 
eralized  residuals)  in  the  MMAE.  Specifically,  as  emphasized  throughout  this  dissertation, 
any  generalized  residual  except  the  traditional  residual  (i.e.,  7  =  1)  should  not  be  used  in 
the  GRMMAE  except  in  the  rare  case  when  the  transformation  matrix  Tj  {ti)  is  the  same 
in  all  elemental  filters,  or  potentially  to  compensate  for  /3-dominance  encountered  using  a 
traditional  MMAE. 

5.1.5  Derived  Equivalenee  Between  MMAE  and  GRMMAE  When  Elemental  Filter 
Have  Same  T  {ti).  The  previous  subsection  mentioned  the  case  when  all  elemental  filters 
have  the  same  transformation  matrix  T  {ti).  This  special  case  applies  to  the  GPS  carrier- 
phase  integer  ambiguity  resolution  application  as  well  as  other  applications  in  which  an 
unknown  measurement  bias  is  estimated.  Section  3.7.1  describes  this  special  case,  demon¬ 
strating  that  the  hypothesis  conditional  probability  for  these  applications  is  independent  of 
the  7  design  parameter.  Specifically,  this  explains  why  the  modified  (i.e.,  post-fit  residual) 
MMAE  successfully  resolved  GPS  carrier-phase  integer  ambiguities,  as  described  in  [17,  18], 
without  being  susceptible  to  /3-dominance.  This  contribution  furthers  the  understanding 
of  the  effect  of  generalized  residuals  on  the  GRMMAE. 
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5.1.6  Derived  Formula  Demonstrating  GRMMAE  Suseeptibility  to  (3-dominanee. 
Using  the  GRMMAE  mathematical  foundation,  the  hypothesis  conditional  probability 
formula,  when  the  transformation  matrix  Tj  (tj)  varies  among  the  elemental  filters,  was 
derived.  This  formula  is  given  in  Equation  (3.40),  repeated  here  for  convenience 


Pki^i)  = 


exp{g*  (U)} 


{Qliti)}  +  Eabsl  g5|jjj||5^exp{(7*(U)} 


(3.40) 


exp 


3  =  1 


where 


Qj  (ti)  =  -  2  K  (^*  )  +  H  (ti)  (Sfc  {ti  )  -  Xj  (t  j)]  1  (f.  )  [rfc  (f.  ) 

+  H  (ti)  (xk  {t- )  -  Xj  {t- ))  1  V  j 


Equation  (3.40)  implicitly  shows  that  the  7  design  parameter,  which  is  part  of  the  trans¬ 
formation  matrices  (U)  and  Tj  (U),  affects  the  (tj)  terms,  but  not  the  exponential 
terms.  This  is  a  direct  cause  of  the  GRMMAE’s  increased  vulnerability  to  /3-dominance  if 
7  =  1  (i.e.,  the  traditional  residual  vector)  is  not  used.  As  described  in  Section  5.1.4,  this 
conclusion  is  an  important  contribution  of  this  dissertation. 

5.1.7  Demonstrated  How 'y  ^  1  Causes  (3- dominanee.  Section  4.3  demonstrated 
in  detail  how  using  7  7^  1  causes  the  GRMMAE  to  be  susceptible  to  /3-dominance.  It 
demonstrated  that  an  improper  choice  for  the  7  parameter  causes  the  /3j  {ti)  terms  to 
become  large  and  increases  the  Euclidean  distance  between  f3*  (ti)  and  /3^  (ti)  for  j  /  k. 
The  conclusion  from  this  section  is  that,  in  general,  7  7^  1  should  not  be  used  in  the 
GRMMAE  as  it  is  currently  implemented.  The  exception  to  this  conclusion  occurs  when 
the  transformation  matrix  T  {ti)  does  not  vary  among  the  elemental  filter.  This  is  a  key 
finding  of  this  dissertation,  as  it  demonstrates  that  the  practice  of  using  post-fit  residuals 
as  in  [17,  18]  is  not  generally  sound. 
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5.2  Recommendation  for  Future  Research 

Due  to  the  GRMMAE’s  susceptibility  to  /^-dominance,  the  GRMMAE  is  not  rec¬ 
ommended  for  widespread  application.  However,  there  are  certain  applications  in  which 
/3-dominance  occurs  with  a  traditional  MMAE.  Most  notably,  /3-dominance  has  been  ob¬ 
served  and  documented  in  the  failure  detection  application  [32,  45].  The  standard  method 
to  address  this  problem  is  to  use  /3  stripping  as  described  in  Section  2.4.1. 

While  /3  stripping  is  an  effective  method  to  overcome  /3-dominance,  it  is  ad  hoc. 
There  is  no  claim  or  guarantee  that  /3  stripping  is  optimal.  However,  the  GRMMAE 
provides  a  mathematical  formulation  in  which  an  appropriate  optimality  criterion  can  be 
defined,  and  a  7  value  chosen  to  achieve  optimal  performance.  Therefore,  a  research  effort 
into  using  the  GRMMAE  in  place  of  /3  stripping  to  address  /3-dominance  is  recommended 
as  a  follow-on  to  this  dissertation.  This  follow-on  research  would  use  the  design  parameter 
7  to  tune  the  (5*  (ti)  terms  to  be  in  consonance  with  the  exponential  terms.  Such  a  study 
is  envisioned  as  follows: 

•  Mathematically  define  the  precise  meaning  for  the  /3t  {ti)  terms  “to  be  in  consonance 
with”  the  exponential  terms.  This  definition  will  likely  quantify  the  relationship 
between  the  magnitudes  of  the  scalar  /3t  {U)  and  exponential  terms  in  the  hypothesis 
conditional  probability  calculation.  It  will  then  probably  set  a  threshold  on  the 
relative  magnitudes  such  that  the  GRMMAE  hypothesis  conditional  probabilities 
are  not  dominated  by  the  /3*  (ti)  terms. 

•  Define  an  optimality  condition,  with  7  as  the  optimization  variable,  to  drive  the 
/3t  {ti)  terms  into  consonance  with  the  exponential  terms. 

•  Derive  an  expression  to  solve  for  the  optimal  7  value.  The  optimal  7  could  be  calcu¬ 
lated  at  each  measurement  update,  or  one  steady  state  value  could  be  determined. 

•  Einally,  the  performance  of  the  GRMMAE  for  /3-dominance  compensation  should  be 
compared  with  the  performance  of  the  traditional  MMAE  implementing  /3  stripping 
to  determine  the  performance  gain  attained  using  the  optimal  method. 

Since  the  7  design  parameter  affects  only  the  /3*  {ti)  term,  there  is  potential  it  can  be 
used  to  compensate  for  /3-dominance.  However,  it  is  not  clear  that  the  use  of  7  will 
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enable  full  “controllability”  of  the  /3-terms  on  all  elementals  filters.  Thus,  implementing 
the  suggested  procedure  is  not  a  trivial  effort.  However,  the  suggested  procedure  is  likely 
at  least  to  provide  an  alternative  to  (5  stripping.  Additionally,  there  is  potential  that  it 
could  provide  improved  performance  and  replace  an  ad  hoc  technique  with  an  optimal  one. 
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